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1 

CHAPTER  1 
INTRODUCTION 

1. 1  Previous  Investigation 

In  this  investigation,  an  analysis  of  the  dynamics  of  bubble¬ 
ring  cavitation  is  performed.  The  analysis  employs  perturbation 
methods  for  analytically  solving  the  governing  nonlinear  ordinary 
differential  equation  for  the  flow.  In  an  effort  to  describe  the 
physical  aspects  of  bubble-ring  cavitation,  the  main  phases  of 
bubble  growth  in  a  region  of  laminar  separation  are  discussed. 

For  a  flow  about  a  hemispherical  headform  containing  a  laminar 
separation  region,  we  assume  the  fluid  contains  a  distribution  of 
nuclei,  invisible  to  the  unaided  eye,  containing  air  and/or  water 
vapor.  These  nuclei  translate  downstream  at  some  velocity  close  to 
the  free  stream  velocity,  V0.  Some  of  these  will  come  in  contact 
with  the  body.  When  a  nucleus  moves  into  the  boundary  layer  on  the 
headform,  it  will  encounter  a  low  pressure  region  which  is  favorable 
to  vaporous  growth.  This  low  pressure  region  has  a  local  static 
pressure  less  than  the  vapor  pressure  within  the  bubble.  The  bubble 
which  had  an  initial  radius,  Rq,  in  the  free  stream  will  grow  to  a 
maximum  radius,  R^,  after  reaching  a  point  on  the  body  where  the 
local  static  pressure  first  equals  the  vapor  pressure.  The  fluid 
conveys  the  bubble  through  the  favorable  pressure  zone  so  that  the 
collapse  phase  commences  and  it  is  here  that  the  maximum  radius 
occurs.  If  there  is  no  separation  on  the  body,  the  bubble  will 
continue  to  collapse  rapidly  and  violently.  For  some  bodies, 
however,  there  can  be  laminar  separation  for  sufficiently  low 
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Reynolds  numbers.  Then  it  seems  possible  Chat  collapse  may  not 
occur  and  the  bubble  will  come  to  rest  within  the  laminar  separation 
region  where  it  undergoes  further  growth  by  diffusion  of  air  from 
the  liquid  into  the  bubble.  This  growth  continues  until  the  bubble 
has  grown  large  enough  to  interact  with  the  free  shear  layer  at  the 
edge  of  the  separation  zone.  This  interaction  causes  the  bubble  to 
translate  downstream  to  the  turbulent  reattachment  region  where  the 
intense  shear  of  this  region  breaks  the  bubble  into  froth.  As  a 
result,  this  froth  creates  a  narrow  ring  of  visible  cavitation  at 
the  downstream  end  of  the  separation  bubble.  This  ring  is  known  as 
bubble-ring  cavitation  which  is  a  form  of  attached  cavitation  and  is 
controlled  primarily  by  laminar  boundary  layer  separation.  Figure  1 
shows  the  configuration  of  a  hemispherical  headform  with  a 
cylindrical  afterbody  and  the  laminar  separation  region  located  on 
the  headform. 

Previous  experiments  have  shown  the  occurrence  of  bubble-ring 
cavitation  to  be  related  to  several  factors.  Holl  and  Carroll  [1] 
observed  the  variation  of  the  influence  of  laminar  boundary 
separation  for  different  test  models  as  a  principal  cause  for 
bubble-ring  cavitation.  For  a  model  whose  configuration  creates  a 
large  separated  region  of  nearly  constant  pressure  with  strong 
pressure  fluctuations  in  the  turbulent  reattachment  region, 
bubble-ring  cavitation  was  observed  at  higher  cavitation  numbers. 
Lowering  the  cavitation  number  brought  on  the  formation  of  a  more 
developed  state  of  cavitation  called  band  cavitation.  Band 
cavitation  is  actually  a  cavity  flow  which  occurs  when  the  laminar 
separation  region  becomes  filled  with  small  attached  cavities. 


Because  bubble-ring  cavitation  serves  as  a  nuclei  source  for  band 
cavitation  it  is  important  to  look  at  the  models  on  which  band 
cavitation  is  also  formed..  Parkin  and  Holl  [3]  observed  band 
cavitation  on  the  hemispherical  nose  and  the  1.5  caliber  ogive  nose. 
Band  cavitation  was  also  observed  on  a  2.0-inch,  1.0  caliber  ogive 
nose  by  Carroll  [2],  a  1.755-inch  diameter  ITTC  nose  by  Arakeri  [4], 
a  2.0-inch,  1/8  caliber  ogive  nose*  by  Keller  [5]  and  a  2.0-inch 
pointed  headfora  used  by  Brockett  [6] .  The  only  models  on  which 
bubble-ring  cavitation  was  observed  were  the  hemispherical  nose  and 
the  1/8  caliber  ogive  nosel  at  various  flow  conditions.  For  a  model 
whose  configuration  has  a  thin  separation  region  with  only  a  slight 
adverse  pressure  gradient  and  small  pressure  fluctuations  at 
reattachment,  bubble-ring  cavitation  did  not  exist  at  all. 

Other  factors  influencing  the  occurrence  of  bubble-ring 
cavitation  are  the  air  content  and  temperature  of  the  water. 

Carroll  [2]  observed  no  bubble-ring  cavitation  for  air  contents 
of  less  than  4.0  ppm.  When  the  air  content  was  held  constant 
at  8.0  ppm,  bubble-ring  cavitation  was  observed  to  disappear 
and  the  limited  cavitation  number,  K£,  decreased  when  the 
temperature  was  increased  at  low  velocities.  Since  raising  the 
water  temperature  decreases  the  solubility  of  air  in  water,  the 
number  and  size  of  the  nuclei  should  also  decrease  suggesting 
that  Kfc  decreases  for  bubble-ring  cavitation.  Carroll  \2] 
observed  this  trend  for  which  suggests  that  air  content  is 
one  of  the  important  factors  controlling  this  phenomenon. 

1.  First  observed  by  Robertson,  McGinley  and  Holl  [20]. 


It  should  be  noted  that  raising  the  temperature  increases  the 
Reynolds  number  so  that  in  marginal  cases  laminar  separation 
is  lost. 

Previous  investigations  by  Arakeri  [71,  Arakeri  and  Acosta 

[8]  and  van  der  Meulen  [9]  of  the  laminar  separation  region  have 

pinpointed  the  initial  separation  point  to  be  downstream  of  the 

minimum  pressure  point  C  .  The  location  of  the  separation 

^min 

zone  did  not  vary  with  Reynolds  number.  Arakeri  and  Acosta  [8] 
were  the  first  to  verify  this  with  Schlieren  photographs  of  the 
thermal  boundary  layer  for  velocities  up  to  60  fps.  Arakeri  [7] 
and  Gates  [8]  experimentally  verified  a  dimensional  variation  of 
the  separated  region  for  various  Reynolds  numbers.  The  tests 
were  performed  on  a  2.0-inch  headform  and  showed  a  decrease  in 
size  of  the  separation  region  with  increasing  Reynolds  number. 

This  dependence,  however,  varies  with  the  shape  of  the  headform. 

This  investigation  requires  a  pressure  distribution  in  the 
region  of  the  separation  bubble  and  the  minimum  pressure  point. 

A  mean  pressure  distribution  for  an  average  free  stream  velocity 
of  40  fps  is  used.  The  experimental  values  for  the  pressure 
coefficient  represent  averages  of  the  measurements  made  by 
Carroll  [2]  and  were  plotted  against  the  dimensionless  axial 
length  X/D.  Since  it  was  desired  to  have  the  pressure  coefficient 
data  plotted  against  the  dimensionless  arc  length  parameter  s, 
a  conversion  was  made  between  X/D  and  s  and  is  shown  in  Appendix  B. 
Figure  2  shows  a  plot  of  the  pressure  coefficient  versus  the 
dimensionless  arc  length  parameter.  Figure  3  shows  the  selected 
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region  of  the  pressure  distribution  used  to  analyze  the  forces  that 
act  on  the  bubbles  and  encoutage  the  occurrence  of  vaporous  bubble 
growth.  Each  figure  shows  the  separation  point  located  downstream 
of  the  minimum  pressure  point  at  a  dimensionless  arc  length 
distance  of  s  -  0.747.  The  location  of  this  separation  point 
agreed  with  theoretical  calculations  made  by  Arakeri  and  Acosta  [8] 
and  van  der  Meulen  [9]. 


1 . 2  Scope  of  this  Analysis 

The  governing  equation  which  describes  isothermal  cavitation 
vapor  bubble  growth  or  collapse,  where  the  static  pressure  is  taken 
to  be  a  function  of  time,  is  the  Rayleigh-Plesset  equation.  Written 
in  dimensionless  form,  excluding  a  viscous  term,  the  isothermal 
Rayleigh-Plesset  equation  for  a  spherical  bubble  is 


rik  +  2<S£)2 
,22 
dx 


1 _ ! 

r3  r 


+  F(ex) 


(1.1) 


This  is  the  form  of  the  equation  used  by  Parkin  [11]  and  is 
the  same  form  used  in  this  investigation.  It  is  a  second  order 
equation  requiring  two  initial  conditions  and  it  is  non- 
autonomous  and  nonlinear.  The  forcing  function  term,  F(ex),  is 
a  time  dependent  pressure  force  which  is  responsible  for  driving 
the  bubble  growth.  The  parameter  e  is  a  dimensionless  small 
parameter  of  the  equation  which  allows  us  to  relate  laboratory 
time  or  "real"  time,  t,  to  a  dimensionless  "bubble"  time  x. 

One  can  write  the  small  parameter  as 


(1.2) 
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where 


T  -  D/V 

o 


(1.3) 


The  parameter  T  is  a  long  time  duration  which  describes  the 
laboratory  time  scale  of  the  forcing  function.  The  definition  of 
T  in  Eq.  (1.3)  allows  the  small  parameter  e  to  scale  automatically 
with  the  value  of  the  free  stream  velocity  V0  and  the  body  diameter 
D.  The  surface  tension  coefficient,  a,  has  units  of  N/m;  the 
density,  p,  has  units  of  Kg/m^;  and  the  initial  nucleus  radius,  Rq, 
has  units  of  meters. 

When  considering  the  effect  of  the  forcing  function,  one  must 
be  aware  of  the  characteristic  time  scales  present  in  this  problem. 
The  laboratory  time,  t,  is  characterized  by  a  slow  time  scale  which 
defines  a  dimensionless  time  duration  across  which  the  forcing 
function  produces  an  environment  favorable  for  vaporous  bubble 
growth.  Thus,  the  dimensionless  slow  time  scale  is  of  order  T  and 
is  written 


Cs  "  T 


(1.4) 


The  second  time  scale  present  in  this  problem  characterizes  the 
individual  bubble  oscillations  that  occur  as  a  bubble  passes 
through  a  varying  pressure  fie.d.  This  fast  time  scale  corresponds 
to  the  dimensionless  bubble  time  t  and  is  written 


Cf  -  T 


(1.5) 


Using  the  definition  of  the  small  parameter  e  as  the  ratio  of  the 
laboratory  time  to  bubble  time,  one  can  write 


Using  Eq.  (1.4),  one  can  write 


which  defines  the  scaling  parameter  for  the  forcing  function  in  the 
vaporous  growth  region.  Relative  comparison  of  the  slow  and  fast 
time  scales  shows  tf  and  tg  to  differ  in  magnitude  by  a  factor  of 
10“ ^  so  that  tf  is  a  very  short  time  duration  compared  to  the  time 
scale  ts  which  defines  the  time  duration  that  the  forcing  function 
acts  in  the  region  of  vaporous  growth. 

Parkin  [11]  derived  an  approximate  parabolic  form  of  the 

forcing  function  making  F(ex)  and  the  governing  equation  non- 

autonomous.  In  an  effort  to  simplify  the  problem,  a  suitable 

combination  of  two  step  functions  was  used  in  place  of  the 

parabolic  form  making  F(er)  piecewise  autonomous  and  making  the 

differential  equation  solvable  by  a  suitable  numerical  method. 

The  choice  of  the  initial  conditions  was  based  on  the  assumption 

that  vaporous  growth  began  with  initial  radius  F  and  ft  *  0. 

o  o 

The  autonomous  form  of  the  isothermal  Rayleigh-Plesset  equation 
was  then  solved  through  the  region  of  vaporous  growth  and  to  the 
maximum  radius. 

In  this  analysis,  the  problem  was  solved  using  a  continuous 
parabolic  representation  of  the  forcing  function  which  was  derived 
from  the  experimental  data  of  Carroll  [2].  Since  the  parabolic 
forcing  function  caused  the  governing  differential  equation  to 
be  nonautonomous ,  the  solution  techniques  used  by  Parkin  [11]  to 


solve  the  piecewise  autonomous  Rayleigh-Plesset  equation  were  not 
fully  applicable.  We  know  that  initially  the  bubble  acts  as  a 
flaccid  bubble  responding  instantaneously  to  the  varying  pressure 
field  encountered  on  the  headform.  The  flaccid  bubble  region  which 
is  primarily  used  to  determine  the  initial  conditions  for  the 
dynamical  problem  has  a  valid  solution  only  up  to  the  initial  point 
where  vaporous  growth  begins.  Beyond  the  flaccid  region  the  bubble 
is  influenced  by  the  inertia  of  the  fluid  surrounding  the  bubble 
as  well  as  the  varying  pressure  field.  The  existence  of  fast  and 
slow  time  scales  in  the  region  beyond  the  flaccid  region  requires 
an  expansion  of  the  governing  differential  equation  as  a  function 
of  the  fast  and  slow  time  scales  as  well  as  the  small  parameter  e. 
This  variation  of  the  method  of  multiple  scales  is  called  the  two- 
variable  expansion  procedure  and  uses  a  perturbation  expansion 
based  on  e  to  separate  out  the  different  solutions  defining  the 
bubble  growth.  Application  of  the  two  variable  expansion  procedure 
produced  a  series  of  nonlinear  differential  equations  and  initial 
conditions  for  only  the  e®  approximation.  As  a  result,  certain 
additional  judicious  approximations  were  required  to  put  the  e® 
solution  in  a  form  which  enables  one  to  tackle  the  e*  and 
systems  of  equations. 


1.3  Motivation  for  the  Investigation 


Previous  experimental  investigations  which  analyzed  bubble¬ 
ring  cavitation  on  hemispherical  headforms  have  revealed  many 
important  details  about  the  particular  flow  phenomena  taking  place. 
As  a  result  of  these  investigations,  the  basic  fluid  dynamics  of 
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this  type  of  flow  are  fairly  well  understood.  On  the  analytical 
side,  most  attempts  at  solution  of  the  governing  equations  are  made 
using  canned  computer  subroutines  that  produce  large  amounts  of 
data  showing  a  detailed  history  of  bubble  growth  or  collapse.  The 
difficulty  with  these  mass  compilations  is  revealed  when  one  desires 
the  effect  of  constraining  a  particular  flow  parameter  at  a 
particular  instant.  Analysis  can  provide  the  option  of  looking  at 
different  classes  of  solutions  and  determining  limiting  cases  of  the 
flow  without  suffering  the  expense,  frustration  and  eventual 
inefficiency  of  large  computer  routines.  The  mathematical 
techniques  being  employed  are  proven  and  work  well  with  these 
complex  flow  equations.  Estimates  of  the  form  of  certain  solutions 
can  often  be  predicted  and  the  existence  of  much  literature  on 
perturbation  methods  is  readily  available. 

Since  the  initial  attempt  to  solve  the  bubble-ring  cavitation 
problem  produced  quite  good  agreement  between  inception  data  and 
theory,  despite  drastic  simplifications,  it  seemed  reasonable  to 
assume  that  if  perturbation  methods  are  used  one  will  produce  equal 
to  better  agreement  between  theory  and  experiment.  It  should  be 
noted  that  even  though  this  method  of  solution  is  more  advanced  than 
most  other  forms  of  analysis,  certain  judicious  approximations  are 
required  to  obtain  a  solution.  An  effort  was  made  to  keep  these 
approximations  to  a  minimum  so  the  true  nature  of  the  flow  can  be 
seen  through  the  solution  of  the  governing  equations.  Successful 
use  of  this  form  of  analysis  may  extend  the  theory's  range  of 
applicability  and  encourage  the  use  of  this  form  of  analysis  for 
other  types  of  cavitating  flows. 
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1.4  Objective  of  this  Investigation 


The  main  objective  of  this  investigation  is  to  find  a  solution 
of  the  nonautonomous  form  of  the  isothermal  Rayleigh-Plesset 
equation  by  use  of  appropriate  perturbation  techniques.  To  do  this, 
a  thorough  parametric  formulation  of  the  forcing  function,  flaccid 
bubble  radius  of  growth,  flaccid  bubble  rate  of  growth  of  the  radius 
and  initial  conditions  across  a  suitable  range  of  cavitation  numbers 
is  performed.  It  was  determined,  for  cavitation  numbers  ranging 
from  K  ■  0.600  to  K  ■  0.700,  that  the  initial  point  of  vaporous 
growth  corresponds  to  dimensionless  arc  length  positions  of 
s  -  0.549  and  s  *  0.600  respectively.  The  initial  conditions,  as 
derived  from  the  flaccid  bubble  relations,  are  applied  at  the 
initial  point  of  vaporous  growth.  The  forcing  function,  which  acts 
across  the  region  of  vaporous  growth,  is  derived  from  the  experi¬ 
mental  data  taken  by  Carroll  [2]  and  is  applied  at  the  initial  point 
of  vaporous  growth.  The  effect  of  the  forcing  function  terminates 
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at  the  separation  point  which  corresponds  to  dimensionless  arc 
length  positions  of  s  ■  0.813  and  s  ■  0.731  for  cavitation  numbers 
of  K  ■  0.600  and  K  *  0.700  respectively.  After  completion  of  these 


tasks,  the  governing  differential  equation  was  expanded  using  the 
method  of  multiple  scales.  Of  the  resulting  set  of  differential 
equations,  the  nonlinear  system  is  solved  and  compared  against 


a  Runge-Kutta  solution  of  the  Rayleigh-Plesset  equation. 
Suggestions  are  made  for  further  study. 
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CHAPTER  2 


FORMULATION  OF  THE  PROBLEM 

2. 1  The  Flaccid  Bubble  Problem 

The  initial  description  of  the  problem  stated  that  the  bubble 
grows  as  a  flaccid  balloon,  experiencing  a  state  of  equilibrium  from 
the  stagnation  point  up  to  the  initial  point  of  vaporous  growth.  In 
order  to  formulate  the  expressions  representing  the  growth  of  the 
bubble,  it  is  necessary  to  discuss  a  few  of  the  fundamental 
assumptions  and  relations  used  in  the  flaccid  bubble  problem. 

To  begin,  we  must  define  a  dimensionless  meridional  arc-length 
parameter  along  the  hemispherical  nose  of  the  body  as 


where  S  is  the  dimensional  arc  length  on  the  body  and  D  is  the 
diameter  of  the  body.  On  the  cylindrical  afterbody,  the  arc  length 
is  the  axial  distance,  X/D.  If  we  assume  the  boundary  layer  to  be 
a  vortex  sheet,  we  can  assume  its  overall  translational  velocity 
to  be  one  half  of  the  local  flow  speed  at  the  edge  of  the  boundary 
layer.  Thus,  it  is  approximated  as 

v (s )  -  A  -  c  (s)  .  (2.2) 

£  p 

Assuming  a  flaccid  bubble  nucleus  which  always  stays  in  the 
boundary  layer,  its  convective  speed  is  size  independent  because 
of  the  assumptions  in  Eq.  (2.2)  but  it  does  change  size 
instantaneously  according  to  the  pressure  variations.  Thus,  we 
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where 


dS  D  dx 

d?  ■  2  or  dS  ■  r  • 


from  Eq.  (2.1).  By  substitution  of  Eq.  (2.2)  and  Eq.  (2.4)  into 
Eq.  (2.3)  and  solving  for  ds/dt,  one  gets 


fr-r^y^  . 


If  we  call  the  parameter  t,  the  laboratory  time,  then  we  relate 
this  actual  time  to  a  dimensionless  bubble  time  by  the  following 
relation 


where  Rq  is  the  nucleus  radius  measured  in  the  free  stream.  The 
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Substituting  Eq.  (2.5)  into  Eq.  (2.7)  and  integrating  over  the 
range  of  experimental  points,  one  gets  the  following  equation 
for  the  dimensionless  bubble  time  along  the  arc  of  the  body 
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The  integral  of  Eq.  (2.8)  depends  on  the  experimental  data  of 
Holl  and  Carroll  (1).  Development  of  a  computer  code  that  could 
produce  an  accurate  value  of  the  integral  would  allow  us  to 
correlate  the  arc  length  parameter  with  the  bubble  time  parameter. 
Care  must  be  taken  when  evaluating  the  integral  numerically  at  the 
stagnation  point  because  it  has  an  integrable  square  root 
singularity.  Evaluation  of  this  singularity  is  discussed  in 
Appendix  A.  The  necessary  requirements  of  a  computer  code  to 
integrate  across  the  range  of  experimental  data  are  the  following: 

(1)  Produce  an  accurate  curve  fit  of  the  integrand 
using  parabolic  interpolation  between  data 
points, 

2 

f(s)  -  As  +Bs+C  . 

where  the  coefficients  A,  B  and  C  are  computed 
for  groups  of  three  successive  data  points. 

(2)  Be  able  to  handle  an  odd  or  even  number  of 
data  points. 

(3)  Integrate  the  parabolic  curve  fit  using  the 
coefficients  A,  B  and  C  for  each  successive 
triad  of  data  points. 

2.1.1  Calculation  of  the  Flaccid  Buddie  Radius,  r(s) 

If  we  consider  the  isothermal  transition  of  the  bubble  from 
the  free  stream  radius  Rq,  to  the  radius  on  the  body  Rj,  then 
we  can  write 


VI  Vi  WWf  HTW  VI  M/l  *C*i  tkr\  nn  m  n  m  n  ■  n 


17 

P,h!  -  PR?  ,  (2.9) 

a  o  g  1 

where  Pa  is  the  partial  pressure  in  the  free-stream  nucleus  and  Pg 
is  the  air  partial  pressure  in  the  bubble  when  it  has  radius  Rj. 
Figure  4  shows  a  typical  nucleus  and  the  internal  and  external 
pressures  that  act  on  the  nucleus.  The  balance  of  pressures  that 
act  on  the  bubble  in  the  free  stream  (R  *  Rq)  is  written  as 

P  +  p  .  |£  +  p  ,  (2.10) 

a  v  R  o 
o 

where  o  is  a  coefficient  of  surface  tension,  Pa  is  the  partial 
pressure  of  the  dissolved  air  in  the  free  stream  nucleus,  P0  is 
the  free  stream  static  pressure  which  establishes  the  external 
environment  for  the  bubble  and  Pv  is  the  vapor  pressure  inside 
the  cavitation  bubble.  On  the  body  (R  »  Rj)  one  can  write  the 
balance  of  pressures  as 

Pg  *  pv  ■  ff+  P<S)  '  (2-U> 

where  Pg  is  the  partial  pressure  of  the  same  mass  of  gas  as  was 
in  the  nucleus  at  R  *  Rq,  but  measured  at  R  =  R^  with  the  size 
of  Rj(s)  being  dependent  upon  its  location  on  the  headform.  The 
external  pressure  term,  P(s),  varies  with  the  arc  length  of  the 
body.  Introducing  a  dimensionless  radius  as 

R1 

r  *  ,  (2.12) 

o 

Eq.  (2.11)  can  be  written  as 


Substituting  Eqs.  (2.10),  (2.12)  and  (2.13)  into  Eq.  (2.9)  one  can 
write 
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(I5-  +  P  -  P  )R3  -  (f^-  +  P(s)  -  P  )r 3R3  . 

o  v'  o  '•R  r  vJ  o 

o  o 


(2.14) 


Using  the  following  definitions  for  the  cavitation  number  and  the 
pressure  coefficient,  one  has 


K  - 


P  -  P 
o  v 


1/2PV 


(2.15) 


P(s)  -  P 


1/2PV' 


(2.16) 


Equation  (2.14)  can  be  written  as 


r+  K(1/2PVo)  -  [^7+  (Cp  +  K)l/2pV^]r 


r2o  1 


2,  3 


(2.17) 


Multiplying  both  sides  of  Eq .  (2.17)  by  Rq/20  and  using  the 
definition  of  the  Weber  number, 


pV2R 

We(Ro’Vo)  * 


(2.18) 


we  write  the  cubic  flaccid-bubble  equation  as 


(%  +  K^We  3  2  KWe 

- r3  +r2  .  i  -™l-0 


(2.19) 


Equation  (2.19)  could  be  solved  exactly,  but  because  r  is  close  to 
unity,  an  approximation  of  the  form, 


rviri*.rv ww ra a#*. wawwivwwwv  •"  *wwnjr  n.*pvw rur FiWnj*i*¥*¥^**ir*rnir*r-nv-*.in 
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r  -  1  -  x  ,  (2.20) 

is  used  where  0  <  x  «  1.  If  one  substitutes  Eq .  (2.20)  into 
(2.19),  he  gets 


A(- 


3xZ  -  3x)  +  x2  -  2x  + 


(2.21) 


where 


(C  +  K)We 

A  -  ■  P-,  — 


(2.22) 


KWe  3 

Because  A  <  — —  and  x  «  1,  the  term  Ax  can  be  neglected  and  the 


remaining  quadratic  equation. 


(3A  +  l)x2  -  (3A  +  2)x  + 


has  a  solution  of  the  form 


C  We 


0  , 


(2.23) 


3A  +  2 
2(3A  +  1) 


[  1  * 


i  -i*L±_ILCWe  j  . 


( 3A  +  2)2  P 


(2.24) 


In  the  limit  as  A  -*■  0,  the  negative  square  root  must  be  chosen 
to  satisfy  the  inequality  0  <  x  <<  1,  therefore 


lim  x  ■  1 
A  -*■  0 


Vi  ^  KWe 
1  +—  • 


(2.25) 


Using  a  binomial  expansion  on  the  square  root  term  in  Eq .  (2.25), 
one  can  write 


wxTrwmramMiwt  nn  vw 
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-V77^ 


■  _  fl  ,  KWe  (KWe)  1 

1  l1  +  ft  ftMM  +  ***J  » 


8  8(16) 


(2.26) 


which  simplifies  the  expression  for  x  to 


\we  /  ,  ,  rave  > 

X  "  T  (“  1  +  16"^  ‘ 


(2.27) 


In  the  limit  as  Cp  ■  -  K  +  0,  Eq.  (2.25)  gives  x  »  0  which 
corresponds  to  a  nucleus  traveling  in  the  free  stream.  The  same 
result  can  be  obtained  directly  from  the  cubic  flaccid  equation 
for  (1)  Cp  -  -  K  -  0  and  (2)  Cp  -  -  K  *  0. 

For  Cp  -  -  K  *  0,  Eq.  (2.19)  simplifies  to 


r2-l 


(2.28) 


which  has  a  root  of  r  »  1  as  we  would  expect.  For  Cp  -  -  K  t  0, 
the  cubic  term  is  neglected  as  before  and  the  equation  for  r  is 
written  as 


2  .  KWe 

r  - 1  -  — 


(2.29) 


or 


V7^ 


1  -  X 


(2.30) 


Therefore, 


/.  .  KWe 

-  v  ~ 


(2.31) 


which  exactly  agrees  with  Eq .  (2.25).  A  check  of  the  validity  of 


Eq.  (2.31)  was  made  by  solving  Eq .  (2.19)  for  and  substituting 


computed  values  of  r  from  Eq .  (2.29)  at  different  cavitation 
numbers.  The  computed  values  of  Cp  were  then  compared  with  the 


iwwrvTv^ 
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experimental  Cp  data.  It  is  clear  from  Fig.  5  that  there  is 
good  agreement  between  the  values  of  Cp  calculated  from  the  cubic 
flaccid  equation  and  the  experimental  data.  Thus,  Eq .  (2.31)  is 
accurate  enough  to  be  a  useful  expression  for  x. 

Calculation  of  r(s)  can  now  be  made  from  the  stagnation  point 
up  to  the  initial  point  of  the  vaporous  growth  region  from 
Eqs.  (2.20)  and  (2.24).  Substitution  of  (2.22)  and  (2.24)  into 
(2.20)  gives  an  expression  for  r  in  terms  of  Cp,  K  and  We.  The 
equation 


r  “  8  +  6(C  +  KjWe  t  1  ±  V  1 


(3(Cp  +  K)We  +  8) 2  P 


(2.32) 

is  plotted  for  various  values  of  K  in  Fig.  6.  Calculation  of  r 
for  various  cavitation  numbers  is  necessary  in  the  derivation 
of  the  initial  conditions  which  is  performed  in  a  later  section. 


2.1.2  Calculation  of  the  Flaccid  Bubble  Growth  Rate,  r(s) 

Another  important  facet  of  the  flaccid  bubble  problem  is  the 
bubble  growth  rate  r(s)  where 


r(s) 


dr  (s ) 


d  x 


(2.33) 


Because  the  bubble  is  conveyed  toward  a  region  where  the  changing 
external  pressure  directly  affects  the  radius  and  therefore  the 
rate  of  growth  of  the  radius,  we  say  the  growth  rate  is  given  by 


j  dC  . 

•  dr  p  ds 

r  dC  ds  dx 

P 


(2.34) 


EXPERIMENTAL  C  DATA 
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DIMENSIONLESS  ARC  LENGTH,  s 

Figure  5.  Comparison  of  the  Experimental  Cp  Data  against  Values 
of  Cp  Calculated  from  the  Cubic-Flaccid  Bubble 
Equation . 


REGION  OF 
INTEREST  FOR 


Figure  6.  Plot  of  the  Dimensionless  Bubble  Radius  versus  the 
Dimensionless  Arc  Length  within  the  Flaccid  Bubble 
Region  for  Cavitation  Numbers  of  K  =  0.60  -*•  K  =  0. 
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Derivation  of  each  component  derivative  in  Eq.  (2.34)  results  in  an 


expression  for  r(s).  Writing  the  cubic  flaccid  bubble  equation  as 


(C  +  K)We 


3  ,  2  .  KWe 

r  +  r  -  1 - 1— 

4 


(2.35) 


and  using  implicit  differentiation  of  Eq.  (2.35),  one  gets 


(C  +  K)We 


,  2  dr  ,  3  We  _  dr 

3r  d^  +  r  r  +  2r  d — 

P  P 


(2.36) 


Solving  for  dr/dCp,  one  has 


' p  3r  (Cp  +  K)We  +  8r 


(2.37) 


which  is  dependent  on  the  experimental  data  and  the  bubble  radius 


r  (s ). 


The  term  dCp/ds  depends  on  the  approximated  parabolic  form  of 


the  experimental  data  which  is 


C(s)*As  +Bs+C 
P 


(2.38) 


If  one  differentiates  Cp(s),  he  finds  that 


*  2As  +  B 


(2.33) 


Finally,  we  know  the  differential  form  of  the  expression  for  the 


dimensionless  time  parameter  is  written  as 


!WW  onwmjrmrTO?  • 
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dt 


ds 


(2.39) 


pR~  'o  /I  -  C  (s) 
K  o  p 


where  the  term  D/V0  -  T  is  the  characteristic  measure  of  laboratory 
time  t,  as  defined  by  Eq.  (1.3).  Referring  to  Eqs.  (1.2)  and  (1.3) 
and  rearranging  for  ds/di,  one  gets 


ds  1 


dx  T 


2o  "  -Cp(s) 


(2.40) 


-  e  /l  -  C  (s) 
P 


where  e  is  defined  by  Eq.  (1.2).  By  substitution  of  Eqs.  (2.37), 
(2.38)  and  (2.40)  into  (2.34),  one  gets 


r(s) 


(2As  +  B)r2We3/2  *^7/2(l  -  C  (s)) 

8  +  3r(C  +  K)We 


(2.41) 


The  sign  of  r(s)  depends  on  what  the  sign  and  magnitude  of 
the  parabolic  coefficients  are.  Calculations  for  r(s)  using  the 
values  of  r(s)  and  the  parabolic  coefficients  up  to  the  beginning 
of  the  vaporous  growth  region  give  the  curves  shown  in  Fig.  7 
at  K*  0.60  K  «  0.70. 
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2.2  The  Forcing  Function 

2.2.1  Definition  of  the  Forcing  Function 

As  we  discuss  the  translation  of  the  bubbles  over  a  headform, 
we  must  be  aware  of  the  driving  force  which  creates  the  environment 
for  subsequent  vaporous  bubble  growth  and  collapse.  We  have 
discussed  the  fact  that  for  vaporous  growth  to  occur  there  must  be 
a  region  on  the  headform  where  the  static  pressure  of  the  flow  is 
less  than  the  vapor  pressure  of  the  water.  This  region  starts  on 
the  headform  at  various  arc  length  positions  depending  on  the 
choice  of  the  cavitation  number.  Downstream  of  the  initial  point 
of  this  growth  region  the  pressure  is  characterized  by  the 
inequality, 

Pv  >  P  >  Pmin  (2*42) 

or 

K  >  -  C  (s)  >  -  C  .  (2.43) 

^  ^min 

Equation  (2.43)  suggests  that  the  negative  of  the  pressure 
coefficient  can  be  used  to  measure  the  force  on  the  bubble  that 
causes  vaporous  growth. 

If  the  bubble  collapses,  then  we  write  the  inequality  as 

K  <  -  C  (s)  .  (2.44) 

P 

Thus,  we  can  set  the  zero  of  the  forcing  function  at 

C  =  -  K  (2.45) 

P 

and  define  the  forcing  function  which  acts  on  the  bubbles  as 


| 


F(er)  where 
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F(e-r)  -  -  C  (b)  -  K  -  -  (C  (s)  +  K)  .  (2.46) 

P  p 

It  should  be  noted  in  this  particular  case  that  the  region  of 
vaporous  growth  and  collapse  extends  along  the  arc  of  the  body  only 
as  far  as  the  separation  point  of  the  laminar  separation  region. 

Once  the  bubbles  enter  the  separation  region  we  assume  that  the 
process  of  vaporous  growth  terminates  and  the  growth  that  ensues 
is  due  to  diffusion  of  air  from  the  liquid  into  the  bubble.  The 
forcing  function  which  is  derived  from  the  pressure  distribution 
can  be  defined  along  the  entire  arc  of  the  body,  but  will  be  used 
in  this  study  only  up  to  the  separation  point.  Figure  8  is  a 
schematic  diagram  of  the  forcing  function  with  its  vaporous  growth 
and  collapse  regions.  The  beginning  of  the  positive  growth  region 
is  located  at  the  intersection  of  the  line  of  Cp  -  -  K  and  the 
forcing  function.  This  initial  point,  where  F  -  0,  is  designated 
as  the  origin  of  the  forcing  function.  The  horizontal  axis  which 
coincides  with  the  line  Cp  -  -  K  defines  the  shifted  er  axis  along 
which  the  duration  of  the  vaporous  growth  process  is  measured.  As 
shown  in  Appendix  B,  this  dimensionless  time  parameter,  which  is 
used  to  scale  the  forcing  function,  is  derived  from  the 
dimensionless  arc  length  parameter  and  has  a  zero  value  at  the 
origin  of  the  er  -  F  coordinates.  Growth  or  collapse  is  designated 
by  et  >  0  up  to  the  point  of  separation  beyond  which  no  more 
vaporous  growth  takes  place. 

When  we  consider  different  flow  conditions  causing  the 
cavitation  number  to  change,  the  line  where  Cp  =  -  K  shifts  up 
or  down  depending  on  the  magnitude  of  K.  For  increasing  values 
of  K,  the  line  of  Cp  =  -  g  shifts  upwards,  decreasing  the  positive 
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vaporous  growth  region  and  increasing  the  negative  collapse  region. 


The  opposite  trend  occurs  when  the  value  of  K  is  decreased  resulting 


in  a  downward  shift  of  the  line  Cp  ■  -  K.  In  view  of  the  location 


of  the  separation  point  and  the  value  of  the  pressure  coefficient 


at  separation,  choosing  K  ■  -  Cp  at  separation  results  in  a 


favorable  region  for  vaporous  growth  which  is  entirely  positive. 


As  stated  earlier,  separation  occurs  at  a  dimensionless  arc  length 


which  corresponds  to  Cp  ■  -  0.6597.  Thus,  positive  growth  exists 
for  all  values  of  Cp  <  -  0.6597  and  growth  and  collapse  occurs  for 


-  0.6597  >  C  >  C  .  Figure  9  shows  the  selected  region  of  the 


P  P, 


forcing  function  calculated  from  the  experimental  data  of  Carroll 


[2]  across  the  region  of  vaporous  growth  for  various  values  of  the 


cavitation  number. 


2.2.2  Axis  Shift  for  the  Forcing  Function 


It  was  seen  that  by  choosing  different  values  of  the 


cavitation  number,  significant  changes  in  the  forcing  function 


take  place  with  respect  to  the  regions  of  growth  and  collapse. 


Accompanying  this  change  in  the  growth  and  collapse  regions  is 


a  shift  in  the  origin  of  the  forcing  function  which  is  defined 


as  F  *  0  at  ex  *  0.  As  the  line  of  Cp  ■  -  K  moves  up  or  down, 


the  origin  of  the  forcing  function  translates  along  both  the  F 


axis  and  the  et  axis.  Figure  8  shows  the  line  along  which  the 


axis  shift  takes  place.  It  is  important  to  know  how  much  the 


origin  translates  for  various  cavitation  numbers  so  that  the 


analysis  can  be  generalized  for  varying  conditions.  For  the 


range  of  cavitation  numbers  being  investigated,  K  *  0.60  to 
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K  ■  0.70,  the  curve  along  which  the  axis  shift  occurs  is  a  straight 
line.  To  define  accurately  the  axis  shift  it  is  necessary  to 
determine  the  variation  of  F  and  er  with  respect  to  the  cavitation 
number . 

From  Figure  10,  the  linear  portion  of  the  forcing  function 
between  the  origin  and  the  line  Cp  ■  -  K  -  -  0.70  is  the  line 
along  which  the  axis  shifts.  The  equation  defining  this  line 
is 

F  -  2.5954( er)  (2.47) 

which  can  be  expressed  in  terms  of  the  cavitation  number  by  the 
relation 

er  -  0.3927600  -  0.23566  .  (2.48) 

Equation  (2.48)  defines  the  er  parameter  along  the  line  of  axis 
shift  for  cavitation  numbers  between  K  ■  0.60  and  K  -  0.70. 
Therefore,  one  can  write 

F  -  1 .01937 ( K)  -  0.61162  (2.49) 

for  K  *  0.60  to  K  *  0.70.  The  reader  may  refer  to  Appendix  B  for 
the  procedure  used  to  calculate  Eqs.  (2.47)  and  (2.48). 

Equations  (2.48)  and  (2.49)  define  the  shift  of  the  origin  for 
different  values  of  the  cavitation  number.  Use  of  the  axis  shift 
is  incorporated  into  the  final  formulation  of  the  forcing  function 
which  is  approximated  by  two  different  parabolic  curve  fits.  It 
should  be  noted  that  the  parameter  ei  is  being  used  because  of 
the  forcing  function  dependence  on  the  single  power  of  the  small 
parameter  e.  Figure  11  shows  a  plot  of  the  forcing  function  for 
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various  values  of  the  cavitation  number  after  the  axis  shift  has 

been  applied*  It  is  clear  from  the  curves  shown  in  Fig.  9  that 

the  axis  shift  places  the  initial  point  of  each  curve  at  a  common 

origin.  Thus,  we  are  able  to  see  clearly  the  effect  of  the  forcing 

function  for  different  cavitation  numbers.  Larger  values  of  the 

cavitation  number  will  cause  the  forcing  function  to  have  a  negative 

effect  near  the  separation  point  while  smaller  values  of  the 

cavitation  number  will  make  the  forcing  function  entirely  positive 

across  the  vaporous  region.  Choosing  a  cavitation  number  greater 

than  or  equal  to  the  negative  of  C  would  result  in  zero 

15  min 

vaporous  growth  for  a  nucleus  as  it  translates  along  the  body. 

2.2.3  Parabolic  Curve  Fit  of  the  Forcing  Function 

Previous  analysis  of  this  problem  was  performed  by  Parkin  [11] 
in  which  the  forcing  function  term  included  within  the  Rayleigh- 
Plesset  equation  was  approximated  by  a  combination  of  step 
functions.  Using  the  definition  of  the  forcing  function  given  in 
Eq.  (2.46),  three  parabolic  approximations  were  made  to  simplify 
the  solution  of  the  Rayleigh-Plesset  equation. 

2. 2. 3.1  Piecewise  Parabolic  Curve  Fit 

In  the  same  manner  that  the  experimental  data  were 
approximated  by  a  piecewise  parabolic  curve  fit,  the  forcing 
function  was  approximated  by  the  same  method.  Using  Eq.  (2.46), 
values  of  F  were  computed  for  various  cavitation  numbers. 

Parkin  [11]  scaled  the  "bubble"  time  parameter  t  against  the 
laboratory  time  parameter  t  by  a  small  parameter  which  we  call 


e,  where 
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Writing  Eq.  (2.8)  as 


where 

Si  . 

Kb)  -  /  —  d!l 

s.  /r-  c  (s) 

1  p 

one  can  rewrite  Eq.  (2.50)  as 


(2.50) 


(2.51) 


(2.52) 


or 

I(s)  -  er  .  (2.53) 

Recall  that  the  laboratory  time  scale  T  across  which  the  forcing 
function  acts  is 

T  «  —  .  (2.54) 

o 

Since  the  small  parameter  e  is  the  ratio  of  bubble  time  to 
laboratory  time,  its  value  gives  the  scaling  factor  between 
the  two  times.  This  scaling  allows  one  the  freedom  of 
choosing  a  different  diameter  headform  and  at  the  same  time 
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maintain  a  proper  scaling  through  the  small  parameter,  e.  Thus 
the  calculated  value  of  the  integral  equals  the  parameter  ct  and 
is  tabulated  in  Appendix  B  with  the  corresponding  values  of  the 
forcing  function.  The  tabulation  of  F  and  ex  was  then  used  to 
interpolate  using  a  piecewise  parabolic  fit  of  the  forcing  function 
for  consecutive  triads  of  points.  The  curve  fit  was  performed  for 
various  cavitation  numbers  between  K  ■  0.60  and  K  ■  0.70.  This 
form  of  the  curve  fit  was  used  to  inspect  carefully  the  region  of 
vaporous  growth  and  collapse  for  different  values  of  K. 

2. 2. 3. 2  Two-Parabola  Curve  Fit 

In  an  effort  to  simplify  the  formulation  of  the  problem,  the 
forcing  function  was  refit  with  two  parabolas  instead  of  a  piece- 
wise  parabolic  fit.  The  two-parabola  fit  incorporated  Eqs.  (2.48) 
and  (2.49)  which  defined  the  axis  shift  for  different  cavitation 
numbers.  Figure  12  shows  the  key  parameters  and  conditions  used 
to  fit  the  two  parabolas  to  the  forcing  function.  The  derivation 
of  the  first  parabola  was  based  on  the  following  conditions: 

(1)  The  parabola  has  the  form 

Fj  -  Aj(ex)2  +  Bj(ex)  +  C j 

(2)  Fj ( ex  -  0)  -  CL  -  0 

(3)  F  (ex  )  =  F  =  A. (ex  )2  +  B,(ex  )  +  C. 

1 v  m'  m  m'  lv  m'  1 


dFl 

(4)  - r  (ex  1  =  2A.(ex  1  +  B  =  0 

d ( ex)  v  m '  1  ^  m'  1 


J 


Conditions  (2)  and  (3)  simply  define  the  value  of  the  forcing 
function  at  the  origin  and  maximum  point  while  Condition  (4)  states 
that  the  forcing  function  has  zero  slope  at  the  maximum  point.  A 
simultaneous  solution  of  these  equations  using  Conditions  (2),  (3) 
and  (4)  gives  the  following  expressions  for  the  parabolic 
coefficients  of  the  first  parabola: 


(ex  y 


1  ex 


Cj  -  0.0 


(2.54) 


The  conditions  used  to  approximate  the  second  parabola  are: 
(1)  The  parabola  has  the  form 


F2  *  A2(ct)  +  B2(ct)  +  C2 


(2)  F-fex  )  *  F  *  A,(et  V  +  bJet  )  +  C.  , 
2'-  m-'  m  2V  m'  ' 
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(3>  dTerT  '  ZAl(eTm)  +  B2 


(4)  F2(ETf)  =  Ff  =  A2(ETfr  +  B,(ETf)  +  C,  . 


2^  f '  2 
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Satisfaction  of  Conditions  (2),  (3)  and  (4)  for  the  second  parabola 
gives  the  parabolic  coefficients 


(2.55) 


Figure  13  shows  the  two-parabola  fit  plotted  against  the 
actual  forcing  function.  To  check  the  validity  of  the  curve  fit, 
a  comparison  of  the  total  impulsive  effect  of  the  forcing  function 
and  the  two-parabola  fit  was  made.  The  impulse  is  written  as 

Jp  ■  /  F( ei)d (ex)  (2. 5b) 

and  was  calculated  using  the  integration  routine  discussed  earlier. 
Figure  14  shows  the  comparison  of  the  impulses  across  the  range  of 
the  vaporous  growth  region.  The  over-estimation  of  the  curve  fit 
resulted  in  a  relative  error  of  3.0391%. 
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2. 2. 3. 3  Trigonometric  Curve  Fit 

An  alternative  to  the  piecewise  parabolic  fit  and  the  two- 
parabola  fit  of  the  forcing  function  is  a  curve  fit  using 
trigonometric  functions  to  produce  a  smooth  continuous  curve  fit 
across  the  entire  range  of  vaporous  growth.  This  type  of  curve 
fit  eliminates  the  need  to  derive  new  initial  conditions  at  the 
junction  between  the  two  parabolas  discussed  in  the  previous 
section.  One  can  assume  that  in  the  interval 

0  <  ts  <  tf  , 

that 

F(ex)  -  F(t  )  -  FfSin(fit  )  (2.57) 

S  L  S 

where  flta  *  it/ 2  when  t3  ■  tf.  Therefore 

q  .11- 

2  Cf 

and 

t 

F(et)  -  F(tg)  -  FfSin(i^)  .  (2.58) 

If  we  consider  the  forcing  function  to  be  the  cause  of  all  primary 
resonances  in  the  system  as  discussed  by  Nayfeh  and  Mook  [16],  one 
writes  the  forcing  function  as 


Fc  Sinftt  =  eKSinfoj  T  +  oT,} 
f  s  o  o  1' 


(2.59) 
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Thus,  the  primary  effect  will  be  in  the  equation  and  we 
maintain  an  autonomous  system  in  the  equation.  That  is  to  say 
the  system  is  soft,  not  stiff.  This  result  is  consistent  with 
the  outcome  of  the  parabolic  forcing  functions  as  can  be  seen  by 
the  inspection  of  those  results. 

2.3  Formulation  of  the  Initial  Conditions 
2.3.1  Derivation  of  the  Initial  Radius,  r(0) 

In  an  effort  to  solve  the  governing  ordinary  differential 
equation  in  the  vaporous  growth  region,  the  initial  conditions  at 
the  beginning  of  this  region  must  be  derived.  The  initial  point 
was  defined  earlier  as  that  point  on  the  body  where  Cp  *  -  K. 
Also,  the  initial  point  is  designated  by  the  dimensionless  time 
t  -  0.  The  conditions  for  the  radius  and  rate  of  growth  of  the 
radius  are  derived  as  functions  of  t  at  r  *  0  and  as  functions 
of  K  across  the  range  K  -  0.60  •>  K  »  0.70.  Thus,  an  arbitrary 
choice  of  the  cavitation  number  would  produce  a  corresponding 
set  of  initial  conditions  for  r(0)  and  r(0). 

Since  the  initial  conditions  are  functions  of  the  parameter 
t,  a  correlation  is  made  between  t  and  s  because  the  initial 
formulation  of  r  and  r  was  made  against  s.  It  was  shown  earlier 
that  by  choosing  different  values  of  the  cavitation  number  the 
axis  of  the  forcing  function  shifted  because  of  the  translation 
of  the  initial  point  of  the  vaporous  growth  region.  To  find 
the  arc  length  positions  corresponding  to  the  initial  point  of 
the  vaporous  growth  region,  one  must  use  the  parabolic  form  of 


the  pressure  coefficient. 
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C  (s)  -  As2  +  Bs  +  C 
P 


along  with  the  condition  that  the  initial  conditions  are  derived  at 
the  point  Cp  -  -  K  from  the  flaccid  bubble  equations.  The  resulting 
equation  is 


As^  +  Bs  +  C  +  K“  0 

which  can  be  solved  with  the  quadratic  formula  to  give 

where 

4A(C  +  K)  .  , 


(2.60) 


(2.61) 


(2.62) 


Equation  (2.61)  is  used  to  compute  the  arc  length  positions  at 
the  beginning  of  the  vaporous  growth  region  which  define  the  range 
across  which  the  equations  for  r  and  r  are  valid.  Using  the 
computed  values  of  r  and  r  at  the  beginning  of  the  vaporous  growth 
region  gives  the  initial  conditions  for  the  problem  at  various 
values  of  K. 

The  initial  condition  for  r  is  equivalent  to  Eq .  (2.30)  which 
was  derived  at  the  point  Cp  *  -  K.  Use  of  a  binomial  expansion 
transforms  (2.30)  into  the  form 


WWW 
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Figure  6  shows  the  region  of  interest  where  the  initial  conditions 
are  sought.  Figure  15  is  a  blow-up  of  that  region  with  the  line 
representing  the  initial  conditions  drawn  through  the  various  r 
curves  at  different  values  of  K.  Equation  (2.63)  is  linear  which  is 
exactly  the  line  of  the  initial  condition  represented  in  Figure  15. 
Thus,  Figure  16  represents  the  initial  condition  r(0)  for  various 
values  of  the  cavitation  number.  In  the  formulation  of  the  problem, 
a  general  initial  condition  of 


r(0)  ■  1  +  q 


(2.64) 


is  used  where  q,  from  Eq.  (2.63),  is  written  as 


KWe 


q  «-g- 


Equation  (2.64)  is  a  more  precise  form  to  use  while  exposing  the 
difference  in  the  initial  condition  used  by  Parkin  [11]  which 
stated 


r (0)  -  1  . 


2.3.2  Derivation  of  the  Initial  Growth  Rate,  r(0) 

In  the  case  of  the  initial  condition  for  r,  the  following 
conditions  were  substituted  into  (2.41)  to  get  an  expression  for 
the  initial  rate  of  growth: 


C  -  -  K 

P 


r  (0) 


'1  + 


KWe 


(2.65) 


.«>  - ir  I- 1  * 


I 


Mil 


mmnm 


w.  UIV.W.TV.  V^.V.V.V.V.V.VIV.v.'ww.-v^v.vi  in 


50 

The  resulting  equation  Is  written  as 

r*(0)  -  |  (1  +  ~-)We3/2  jj2-  -^B2  -  4A(C  +  K)  (2.66) 

and  Is  a  function  of  the  cavitation  number  and  the  physical 

parameters  of  the  flow  and  headform.  Figure  17  shows  the  line  of 

initial  condition  for  r(0)  across  the  range  of  cavitation  numbers 

being  investigated.  In  bubble  time  the  initial  condition  has  a 

constant  value  with  a  magnitude  of  the  order  1  x  10  \  In 

laboratory  time  this  translates  to  a  velocity  with  a  magnitude  of 
-9 

the  order  1  x  10  fps.  Since  the  magitude  is  relatively  small, 
the  initial  condition  for  r(0)  is  approximated  as  zero,  although 
the  fact  that  it  is  not  actually  zero  is  found  to  be  conceptually 
important  later.  The  zero  value  for  r(0)  corresponds  to  the 
initial  condition  used  by  Parkin  [11]. 

In  summary,  using  the  properties  of  a  flaccid  bubble  under¬ 
going  isothermal  expansion,  expressions  were  derived  for  the 
radius  of  growth  and  the  rate  of  growth  of  the  radius.  These 
expressions  are  functions  of  the  cavitation  number,  the  particular 
parabolic  coefficients  used  to  approximate  the  experimental  data 
and  the  physical  parameters  of  the  flow  and  the  headform.  These 
expressions  are  valid  in  the  flaccid  bubble  region  only,  beginning 
at  the  stagnation  point  and  ending  where  the  pressure  coefficient 
equals  the  negative  of  the  cavitation  number,  i.e.,  Cp  =  -  K. 

These  values  of  r  and  r  at  the  point  C„  -  -  K  represent  the 


5g**rf~le^ 


initial  conditions  used  to  solve  the  governing  differential 
equation.  The  resulting  initial  conditions  are 


r(0)  -  1  +  Q  ,  Q  <  1  [from  Eq.  (2.64)] 
and  (2.67) 

r(0)  ■  0  [evaluation  of  Eq.  (2.66)] 

for  the  range  of  cavitation  numbers  under  investigation. 


2.4  The  Differential  Equation  for  an  Isothermal  Bubble 

Written  in  dimensionless  form,  the  governing  equation  for 
vaporous  growth  and  collapse  of  a  spherical  isothermal  cavitation 
bubble  is 


r 


2  fdr^ 
2 


j-  +  F(et)  . 


(2.68) 


Equation  (2.68)  is  a  second  order  nonlinear  ordinary  differential 
equation  requiring  two  initial  conditions  for  its  complete 
solution.  The  initial  conditions  for  r(0)  and  r(0)  were  derived 
in  the  previous  section,  Eqs.  (2.63)  and  (2.66).  The  forcing 
function,  F(ex),  which  causes  the  bubble  to  grow  from  the  initial 
condition  at  r(0),  to  some  maximum  radius  rm,  was  approximated  by 
two  parabolas  making  Eq .  (2.68)  nonautonomous .  Past  studies  of 
the  nonautonomous  form  of  Eq.  (2.68)  have  been  performed  using 
appropriate  numerical  methods.  In  several  unpublished  numerical 
studies  by  Parkin,  it  was  found  possible  to  distinguish  four 
classes  of  solution  for  Eq .  (2.68).  A  Class  1  solution  is 
characterized  by  small-scale  oscillations  that  make  up  the  major 


part  of  the  bubble's  motion.  In  a  Class  2  solution,  the  bubble 
motion  involves  a  periodic  motion  similar  to  a  Class  1  solution 
but  with  a  larger  amplitude.  A  Class  3  solution  has  periodic 
solutions  for  the  bubble's  oscillations  but  with  larger  amplitudes 
The  size  of  the  amplitudes  was  found  to  be  dependent  on  the  key 
parameters  of  the  flow.  Finally,  a  Class  4  solution  represented  a 
bubble  growing  infinitely  large.  Of  interest  here  are  the  Class  1 
Class  2  and  Class  3  solutions.  It  is  desired  that  a  solution  of 
Eq.  (2.68)  be  obtained  that  would  allow  a  parametric  study  of 
these  different  classes  of  solution  to  be  performed.  With  the 
initial  conditions  and  forcing  function  previously  derived,  one 
can  now  apply  the  method  of  multiple  scales  in  an  attempt  to  get 
an  approximate  analytic  solution  of  the  isothermal  Rayleigh- 
Plesset  equation. 
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CHAPTER  3 

PARTIAL  SOLUTION  OF  THE  DYNAMICAL  PROBLEM 


3. 1  The  Method  of  Multiple  Scales 

Having  established  the  governing  differential  equation  and 
initial  conditions  for  isothermal  cavitation  bubble  growth  and 
collapse,  the  method  of  multiple  scales  is  used  to  find  an 
approximate  analytical  solution  to  the  problem.  The  method  of 
multiple  scales  is  used  because  of  the  two  dissimilar  time  scales 
that  exist  in  this  problem.  Since  there  are  two  different  time 
scales,  a  variation  of  the  method  of  multiple  scales  called  the 
two-variable  expansion  method  is  used.  Application  of  the  two- 
variable  expansion  method  begins  with  the  expansion  of  the  two 
time  scales. 

3.1.1  The  Time  Scales 

As  discussed  earlier,  there  are  two  time  scales  that 
characterize  this  problem.  The  time  measured  in  the  laboratory  is 
the  slow  time  scale  ts  and  is  the  characteristic  time  scale  of  the 
forcing  function.  The  fast  time  scale  tf  is  a  very  short  time 
compared  to  the  slow  time  scale  and  is  the  characteristic  time 
scale  of  the  individual  bubble  oscillations.  Expansion  of  the 
time  scales  based  on  the  small  parameter  e,  defined  by  Eq.  (1.2), 
is  derived  by  Cole  and  Kevorkian  [14],  The  slow  time  scale  is 
simply 

tg  -  ex  (3.1) 


while  the  fast  time  scale  is  written  as 
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Because  the  individual  bubble  oscillations  are  very  fast,  the 
expansion  for  tf  includes  additional  terms  which  allow  for  different 
frequencies  of  oscillation  that  may  occur  as  the  bubb1 e  passes 
through  a  varying  pressure  field.  Based  on  these  two  time  scales, 
the  dynamical  equation  can  be  expanded  to  the  order  e^. 


3.1.2  Formulation  of  the  Dynamical  Equations  to  Order 

Formulation  of  the  dynamical  equations  to  order  is 
accomplished  by  using  the  expansions  for  the  slow  and  fast  time 
scales,  Eqs.  (3.1)  and  (3.2),  to  derive  the  first  and  second 
derivative  expansions  as  functions  of  the  dimensionless  bubble 
time  r.  One  can  write 

T  -fCVbf)  (3.3) 


where  the  derivative  with  respect  to  Eq.  (3.3)  is 


d _ 

d  t 


( 1  +  e  oi. 


3 

+  e  a), 


(3.4) 


It  follows  from  Eq.  (3.4)  that  the  second  derivative  is 


2  3  2  2 
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f  s 


+ 


(3.5) 


or  written  in  ascending  power  of  e  the  second  derivative  expansion 
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(3.9) 


If  one  assumes  Che  bubbles  under  consideration  are  spherical, 


then  the  term  S(r,n)  from  Eq.  (2.68)  is  set  equal  to  unity.  The 
term  F(ct)  from  Eq.  (2.68)  has  the  form  of  a  parabola  as  discussed 
in  Section  (2.2).  Unfortunately,  use  of  a  parabola  of  the  form 

F  ■  Ae^t^  +  Bet  +  C  (3.10) 

s  s 

requires  the  solution  of  nonautonomous  first  and  second  order 
expansion  equations  to  get  a  complete  solution.  By  using  a 
trigonometric  form  of  the  forcing  function  represented  by 
Eq.  (2.59),  one  needs  to  solve  only  the  nonautonomous  first  order 
expansion  equations  since  the  trigonometric  form  is  expressed 
only  to  the  first  order  of  the  small  parameter,  e. 

Calculation  of  each  term  in  Eq.  (2.68)  based  on  the  general 
perturbation  expansion  for  r  and  the  first  and  second  time 
derivative  expansions,  permits  one  to  write  the  isothermal 
Rayleigh-Plesset  equation  in  ascending  powers  of  the  small 
parameter  e  as  a  function  of  the  fast  and  slow  time  variables. 

The  same  procedure  can  be  performed  on  the  initial  conditions. 
Introduction  of  a  normalized  radius  u  where 


allows  one  to  write  the  normalized  isothermal  Rayleigh-Plesset 
equation  with  initial  conditions  consistent  with  those  used  by 
Parkin  (11].  The  resulting  set  of  differential  equations  and 
initial  correlations  up  to  can  be  solved  and  substituted  back 
into  Eq.  (3.7)  using  (3.11)  to  get  the  final  solution  for  the 
dimensionless  bubble  radius  as  a  function  of  time. 
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The  normalized  equations  and  initial  conditions  written  in 
ascending  powers  of  c  are: 

Order 
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du 
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Uq(0)  -  1.0 


du( 

dt. 


Initial  Conditions 
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Initial  Conditions 
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Order 
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u3(0)  =  0.0 


du^(0)  du^O) 

dtf  ~  w2  “dt^ 


dun(0) 


Initial  Conditions 


3  dt. 


(3.14b) 
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3.2  Solution  of  the  Multiple  Scale  Equations 
3.2.1  The  Solution 

The  first  step  in  obtaining  the  complete  solution  of  the 
isothermal  Rayleigh-Plesset  equation  is  to  solve  the  normalized 
zeroth  order  differential  equation  that  results  from  the 
application  of  the  method  of  multiple  scales.  The  zeroth  order 
equation,  Eq.  (3.11a),  is  still  nonlinear  but  is  now  autonomous 
because  the  time  dependent  forcing  function  has  been  reduced  to 
a  constant.  Equation  (3.11a)  is  analogous  to  the  autonomous 
form  of  the  Rayleigh-Plesset  equation  studied  by  Parkin  [11] 
where  the  term  Fg  corresponds  to  his  piecewise  autonomous  step 
function  representation  of  the  forcing  function.  From  Eq.  (3.11a) 
one  can  examine  the  intricacies  of  this  problem  by  investigating 
the  energy  curves  and  phase  plane  trajectories.  This  form  of 
analysis  is  restricted  to  the  autonomous  zeroth  order  differential 
equation  and  is  very  useful  in  determining  the  limits  of  the 
periodic  and  non-periodic  solutions  which  are  based  on  the 
location  and  character  of  the  singularities  that  exist  in  the 
autonomous  system. 


3 . 2 . 1 . 1  The  Potential  Energy  Function 


The  zeroth  order  equation  and  initial  conditions  may  be  written 


2  2 

u0  +  2  lUI)  =1 _ l _ 1 _ 1 _ + _ -£ _ 

dt2  2  dtf  U2  (1  +  Q)5  U0  (1  +  Q)3  (1  +  Q)2 


uQ(0)  =  1.0 


L_ 


li 


11  (.•  <  *  r  t  l4Ji»  *  i  >  U  l».  IVA  l».  lV*llkra<U6l». 


WTV!T7ntnXn^7ntn«cnK  '.xjih.-k.-mi  •.uni. 


(0)  -  0.0  . 


It  is  of  great  interest  to  consider  the  different  solutions  of 


the  zeroth  order  equation  for  the  two  cases  when  (1)  Fc  •  0  and 


(2)  Fq  *  0.  One  can  replace  the  second  order  differential 


equation  by  a  pair  of  coupled  first  order  equations  using  the 


relation 


dt 


To  simplify  the  expressions  by  dropping  the  subscript  zero, 


the  zeroth  order  equation  and  initial  conditions  can  be  rewritten 


dv  .  3  2  y  1 

u  XT-  +  7  v  *  i - 

dtf  2  u3  (1  +  Q) 


l 

5  ~  u 


(1  +  Q)3  (1  +  Q)2 


(3.15) 


u(0)  -  0 


v(0)  -  1  . 


Using  the  transformation 


dv  3  2  Id  (  3  2-> 

u  dF  +  2  v  ■  77  ST  (“  v ) 

2u 


(3. 1'6) 


one  may  express  Eq .  (3.15)  in  integral  form  as 


32  ,  r2v  1 

u  v  „  /  [_L  - - 


2y  1  2u  2U  FC  i 

— L - t - 7  +  - Hdu 

u  (i  +  q)  (i  +  or  (i  +  or 


(3.17) 
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The  left-hand  side  of  Eq .  (3.17)  is  proportional  to  the  kinetic 
energy  of  the  bubble  motion  while  the  right-hand  side  corresponds  to 
the  negative  potential  energy,  -  V,  of  the  bubble.  Evaluating  the 
integral  in  Eq .  (3.17)  as  an  indefinite  integral,  one  can  write 


-  V 


2y£nu 


(1  +  Q)  (1  +  Q) 


.  2  3 

3  +  3  u 


+  k 


(1  +  Q)‘ 


(3.18) 


where  k  is  a  constant  of  integration  that  permits  one  to  adjust  the 
level  of  V  for  different  initial  conditions.  Setting  Fc  ■  0  and 
evaluating  Eq.  (3.18)  at  the  initial  conditions  defined  for  the 
zeroth  order  equation,  one  writes  the  potential  energy  function  V 


li. 


£nu  + 


1 


(1  +  Q)' 


(1  +  Q)' 


(u2  -  1) 


(3.19) 


and  the  plot  of  Eq.  (3.19)  versus  the  normalized  radius  u  is  shown 
in  Fig.  18  for  values  of  Q  ■  0.30  to  Q  ■  0.0.  It  is  easy  to  see 
from  Fig.  18  that  as  the  value  of  0  decreases  toward  zero,  the  point 
of  minimum  potential  energy  shifts  to  the  right.  The  relationship 
defining  the  location  of  this  minimum  energy  is  derived  by  setting 
the  first  derivative  equal  to  zero.  Thus, 


dV 

du 


2y 


I  + 


2u 


(1  +  Q)2  U  (1  +  Q)3 


0 


or 


_/y_ 


1  +  Q 


(3.20) 
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which  defines  analytically  '.he  location  of  minimum  energy.  In  this 
study  we  restrict  our  attention  to  positive  values  of  u  and  thus 
neglect  the  negative  root.  The  importance  of  Eq.  (3.20)  will  become 
evident  later  when  we  study  the  location  and  character  of  the 
singular  points  of  Eq.  (3.15).  Moreover,  the  translation 
experienced  by  the  minimum  energy  point  is  the  critical  condition 
that  allows  one  to  distinguish  between  two  kinds  of  motion 
characteristic  to  the  autonomous  system. 

If  one  evaluates  Eq.  (3.18)  using  the  initial  conditions  from 
Eq.  (3.15)  and  sets  Fg  *  0,  the  resulting  potential  energy  function 
is  written  as 

V  - - ^ — r  £nu  + - - - s-  (u2  -  l)  -  i - - — j  (u3  -  l)  . 

(1  +  Q;  (1  +  Q)  (1  +  Q) 

(3.21) 

Equation  (3.21)  represents  the  potential  energy  function  for  the 
autonomous  system  having  a  non-zero  forcing  function  and  is  plotted 
in  Fig.  19.  Clearly,  the  contribution  to  the  potential  energy 
function  of  the  last  term  in  Eq .  (3.21)  becomes  overwhelming  as  u 
increases,  causing  the  energy  curve  to  turn  downwards.  As  before, 
the  location  of  the  minimum  and  maximum  points  is  determined  by 
setting  the  first  derivative  equal  to  zero.  Performing  the  same 
operation  on  Eq.  (3.21),  we  find  the  result  to  be  a  cubic  equation 
of  the  form 


3 


u 


V1 


+  Q)' 


Fc(l  +  Q) 


0  . 


(3.22) 
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Figure  19.  Potential  Energy  Plot  for  Fc  f  0. 
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The  roots  of  Eq.  (3.22)  are  found  using  Tartaglia’s  Method. 
For  the  case  of  three  real  and  »  nequal  roots,  the  roots  of 
Eq.  (3.22)  are  of  the  form 


P 

Xi  -3 


(3.23a) 


where 


.  2ir(i  -  1)> 
cos (e1  + - ^ - J 


i  -  1,2,3 


(3.23b) 


1  -1  r3b  n 

0.  -  -r  COS  ( - J 

1  3  vam; 


m  -  2 


a  ■  t  (3q  -  p2) 


b  -  (2p3  "  +  27r) 


The  parameters  p,  q  and  r  are  the  coefficients  of  a  cubic  equation 
having  the  form 


3  2 

u  +  pu  +  qu  +  r  -  0 


(3.24a) 


Therefore, 


P  "  "  Fc(l  +  Q)  ’ 


q  -  0 
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and 


(3.24b) 


Thus,  in  terms  of  the  original  paramo of.  the  cubic  coefficients, 
the  roots  of  Eq .  (3.22)  are 


3Fcd 


+  Q) 


3FC(1 


+  C'i 


c  os ,  fl  , 


2m(i  -  1) 


) 


i  -  1,2,3 


(3.25a) 


where 


■j  COB 


- 1 


27F  ly 
[1 - 3s-] 


•  Evaluating  the  roots  for  Q  »  0.3,  y  -  1.4,  Fq  ■  0.2,  one  gets 
i  -  1  ,  -  3.600356643 


i  -  2  ,  u2  -  -  0.825817193 


(3.25b) 


i  -  3  ,  u3  -  1.071614395 

Since  we  are  interested  in  positive  roots  only,  the  second  is 
neglected  leaving  the  first  and  third  root  to  define  the  location 
of  the  maximum  and  minimum  points  in  Fig.  19. 


3.2. 1.2  Singular  Points  and  the  Phase  Plane 

Having  derived  an  expression  for  the  potential  energy  of  the 
autonomous  system,  we  can  now  look  at  the  phase  plane  in  conjunction 
with  the  potential  energy  to  study  the  singularities  that  exist  in 
the  system.  To  begin  the  analysis  we  need  to  apply  Eq.  (3.16)  to 
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(3.15)  and  use  a  generalized  initial  condition  for  the  initial 
radius  that  states  u(0)  ■  u0.  The  initial  condition  of  the  bubble 
wall  velocity  will  remain  zero  because  its  actual  initial  value  is 
so  small  that  all  analyses  whether  global  or  specific  will  be 
accurately  represented  by  u(0)  -  0.  Use  of  the  generalized 
initial  condition  on  the  bubble  wall  radius  will  enable  us  to 
perform  a  global  phase  plane  analysis  to  explore  the  properties 
of  the  autonomous  system  and  to  develop  some  knowledge  about  the 
physical  parameters  that  lead  to  different  classes  of  solutions. 
Afterwards,  the  initial  conditions  defined  by  Eq.  (3.11b)  will  be 
used  to  look  at  a  specific  set  of  trajectories  in  the  phase  plane 
and  how  they  relate  to  the  global  analysis. 

Applying  Eq.  (3.16)  to  (3.15),  one  can  write 


d  f  3  2-I 

sr  (u  v  ) 


2y  1  2u  .  „  2  FC 

— I - - - -  +  2u  - T  . 

U  (1  +  Qr  (1  +  Q)  (1  +  Q) 


(3.26) 


Using  Eq.  (3.26),  one  can  derive  an  expression  for  the  phase  plane 
trajectory  by  computing  the  derivative  on  the  left-hand  side  and 
solving  the  equation  for  dv/du.  The  result  is 


d  v 
du 


u3F, 


(1  +  Q)3  (1  +  Qr  (1  +  Q) 

4 

u  v 


3  2  3 

2  "  2  V  U 


P(u,v) 

Q(u,v) 


(3.27) 


which  defines  the  bubble  wall  velocity  as  a  function  of  the  bubble 
wall  radius  and  the  other  physical  parameters  of  the  flow.  The 
singular  points  of  the  autonomous  system  are  located  at  those 
points  which  satisfy  the  conditions  Q  ■  P  ■  0.  Because  u  >  0, 

Q  *  0  when  the  initial  condition  for  the  bubble  wall  velocity  is 
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satisfied.  Thus,  all  singular  points  are  located  along  the  u-axis 
at  positions  satisfying  P(u,0)  ■  0.  Evaluating  P(u,0)  ■  0,  one 
gets  the  cubic  equation 


3 

u 


1  2 
Fc(l  +  Q)  U 


+ 


X 

Fc(l  +  Q)3 


0 


which  is  exactly  the  same  as  Eq.  (3.22).  As  we  saw  earlier,  the 
roots  of  Eq.  (3.22)  define  the  locations  on  the  u-axis  of  the 
maximum  and  minimum  points  of  the  potential  energy  function.  Thus, 
we  can  say  that  the  location  of  the  maximum  and  minimum  points  on 
the  potential  energy  curve  correspond  exactly  to  the  location  of 
the  singular  points  of  the  system. 

In  order  to  determine  the  character  of  these  singularities 
one  can  apply  Liapunov's  Method  (see  Ref.  [15))  which  requires 
finding  the  characteristic  roots  of  the  Jacobian  matrix 

I?  (u*o)  <u>o) 

The  value  of  u  in  the  argument  of  each  derivative  function 
corresponds  to  the  values  of  the  real  roots  denoted  by  Eq.  (3.20) 
for  ■  0  and  Eq.  (3.25)  for  Fq  *  0. 

When  Fq  *  0,  the  characteristic  equation  is 

X2  -  -  2u5  (3.28a) 


or 


1,2 


±  i 


/2 


(1  +  Q) 


(3.28b) 
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When  Eq.  (3.28b)  is  evaluated  at  the  positive  root  of  Eq.  (3.20), 
the  result  is  two  purely  imaginary  characteristic  roots  meaning 
the  positive  root  is  characterized  as  a  vortex  point.  In  a  small 
neighborhood  in  the  phase  plane,  all  motion  will  circle  around  this 
vortex  point  thus  defining  a  periodic  solution  for  u. 

When  Fc  *  0,  the  characteristic  equation  becomes 

3F 

,2  2  5  ,  JrC  6 

A  - - -  u  +  - -  u  .  (3.29) 

(1  +  Q)5  (1  +  Q; 

By  evaluating  Eq.  (3.29)  at  the  positive  roots  of  Eq.  (3.22),  one 
can  determine  the  nature  of  those  roots.  For  the  smallest  positive 
root  in  Eq.  (3.25b),  the  characteristic  roots  are  purely  imaginary 
meaning  the  root  is  a  vortex  point.  For  the  largest  root  in 
Eq.  (3.25b),  the  characteristic  roots  are  both  real  meaning  the 
root  is  a  saddle  point.  The  location  of  the  saddle  point 
corresponds  to  the  outermost  limit  on  the  u-axis  that  a  trajectory 
representing  a  periodic  solution  will  pass.  Ma  and  Wang  [17] 
obtained  similar  results  for  this  specific  case.  These  results 
are  also  completely  consistent  with  those  of  Parkin  [11]  who  found 
for  Fq  *  0  that  only  one  vortex  point  exists.  And  for  *  0, 
a  vortex  and  saddle  point  exist  on  the  u-axis  with  the  vortex 
point  located  closer  to  the  initial  radius  u0  than  the  saddle 
point.  Moreover,  the  roots  that  lie  on  the  negative  u-axis  are 
not  pertinent  to  our  problem  and  should  not  affect  the  solution. 
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Rearranging  Eq.  (3.17)  for  the  bubble  wall  velocity,  one  can 

write 


which  defines  the  family  of  trajectories  in  the  positive  phase  plane. 

The  family  of  curves  resulting  from  Eq .  (3.30)  is  a  series  of  closed 

loops  and  open  curves  representing  periodic  and  non-periodic  solutions, 

respectively.  The  plus  and  minus  designation  of  Eq.  (3.30)  defines  the 

trajectories  above  and  below  the  u-axis.  The  extent  of  the  closed  loop 

periodic  solutions  is  determined  by  the  location  of  the  singularities 

which  were  found  to  coincide  with  the  position  on  the  u-axis  of  the 

minimum  and  maximum  potential  energies.  Plots  of  the  potential  energy 

and  phase  plane  trajectories  are  shown  in  Figs.  20  and  21,  Figure  20 

corresponds  to  the  potential  energy  and  phase  plane  trajectories  when 

Fc  ■  0,  y  -  1.4.  Notice  the  location  of  the  minimum  potential  energy 

corresponds  to  the  vortex  point  in  the  phase  plane  about  which  the 

periodic  trajectories  are  focused.  It  is  important  to  point  out  that 

the  location  of  the  minimum  potential  energy  translates  from  left  to 

right  as  the  value  of  Q  decreases  from  Q  ■  0.30  to  0  =  0.00.  The 

dotted  curve  in  the  upper  potential  energy  plot  corresponds  to  the 

case  where  the  trajectory  and  the  vortex  are  isolated  at  the  origin 

in  the  phase  plane  for  a  critical  value  of  Q  .  ■  0.  183216.  This 

crit 

critical  motion  is  then  defined  by  u(t)  -  1  and  u(t)  =  0  for  all  time. 
The  importance  of  this  critical  condition  is  that  it  separates  two 
physically  distinct  types  of  motion  that  the  bubble  experiences 
depending  on  its  initial  size.  The  trajectories  to  the  left  of 


NORMALIZED  BUBBLE  WALL 


Figure  21.  Phase  Plane  and  Level  Line  Potential  Energy  Plot 
for  the  Zero-Order  Analysis  when  F  4  0. 
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u(t)  ■  1  represent  small  scale  oscillations  of  a  flaccid  air  bubble 
while  those  to  the  right  of  u(x)  ■  1  have  larger  amplitudes  of 
oscillation  involving  vaporous  growth.  Since  we  are  considering 
primarily  the  vaporous  growth  region,  we  will  limit  our  discussion 
and  investigation  to  those  trajectories  that  lie  to  the  right  of 
u(t)  ■  1.  The  location  of  the  vortex  points,  which  coincides  with 
the  minimum  potential  energies,  is  defined  by 


u 

v 


(1  +  0) 


(3.31) 


Thus,  one  can  write  at  uv  ■  1 


9crlt  '  ^  '  1 


(3.32) 


for  any  value  of  y,  If  Q  -  0,  then  y  *  1  which  is  the  lowest 
possible  value  of  the  dissolved  air  content  for  positive  Q  values. 
Since  Q  was  previously  defined  in  Eq.  (2.63)  as 


the  minimum  limit  on  Q  corresponds  to  V0  *  0  or  K  *  0. 

If  one  now  considers  the  potential  energy  and  phase  plane 
trajectories  for  *  0  at  y  *  1.4,  it  can  be  seen  from 
Fig.  21  that  the  nature  of  the  curve  is  quite  different  from 
the  case  when  Fq  -  0.  To  simplify  the  analysis  we  will  look 
at  only  one  energy  curve  with  its  corresponding  phase  plane 
trajectories.  From  Fig.  21  one  can  see  that  the  minimum  point 
on  the  energy  curve  coincides  with  the  vortex  point  in  the 
phase  plane  while  the  maximum  energy  point  concides  with  the 
saddle  point.  The  location  of  the  saddle  point  defines  the 
maximum  normalized  bubble  radius  through  which  a  trajectory 
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representing  a  periodic  solution  passes.  Therefore,  all 

trajectories  representing  periodic  solutions  lie  to  the  left  of 

the  saddle  point  and  circle  about  the  vortex  point  located  just 

to  the  right  of  some  generalized  initial  condition  u(0)  *  1  +  Q. 

The  outermost  trajectory  representing  a  periodic  solution  and 

passing  through  the  saddle  point  is  called  the  separatrix.  The 

parameters  that  define  the  separatrix  as  the  outermost  limit  of  all 

periodic  trajectories  is  characterized  by  the  critical  parameters 

y  and  F  .  These  critical  parameters  are  determined  by 
C  crit 

simultaneous  solution  of 

V(u)  -  0  (3.33) 

and 

£  (u)  -  0  .  (3.34) 

Equations  (3.33)  and  (3.34)  are  to  be  evaluated  at  the  location  of 
the  saddle  point  as  determined  from  Eq.  (3.25).  All  trajectories 
that  lie  outside  of  the  separatrix  represent  bubbles  which  will 
grow  to  an  infinite  radius  across  an  infinite  time  interval.  In 
this  study  we  are  mainly  interested  in  the  trajectories  that 
represent  periodic  solutions. 

To  solve  for  the  critical  parameters  that  determine  the 
separatrix,  one  can  write  Eqs.  (3.33)  and  (3.34)  as 

V(u) - £nu - - - T  (u2  -  1)  +  | - — r  (u3  -  l)  -  0 

(1  +  Q)  0  +  Q)  (1  +  Q) 


and 


(3.35) 
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dV 

du 


2y  1  u  FC  2  „ 

- ! - - - -  +  - -  U  -  0 

(i  +  o)  u  (i  +  or  (i  +  o) 


(3.36) 


Solving  Eqs 
for  Fc 

crit 


'crit 


.  (3.35)  and  (3.36)  for  y  and  then  solving  simultaneously 
one  gets 


3(u^  -  1 1  -  6u^lnu _ 

(l  +  Q)[2(u^  -  l)  -  6u^£nu] 


(3.37) 


which  defines  the  critical  forcing  parameter  as  a  function  of  0. 

Rearranging  Eq.  (3.36)  for  y  and  substituting  Eq.  (3.37)  for 

F_  one  can  write  the  expression  for  y  .  as 
Ccrit  crlt 

2,1  „  %  1 

,  o  u  (-r  -  £nuj  - 

■  °  + «' «'  (‘  -  «(-5tt - : — r>]  •  <3-38) 

u  ( J  -  £nu)  -  — 


Thus,  Eqs.  (3.37)  and  (3.38)  are  the  critical  parameters  which 

when  used  in  conjunction  with  Eq .  (3.30)  define  the  critical 

trajectory  called  the  separatrix.  Parkin  [11]  asserted  that  the 

separatrix  corresponded  to  a  barrier  between  Class  1  and  Class  4 

solutions  for  a  non-zero  forcing  function  parameter,  Fc.  Recall 

that  a  Class  4  solution  corresponds  to  a  bubble  which  grows  to  an 

infinite  radius  in  an  infinite  amount  of  time.  If  one  chooses  a 

value  of  the  forcing  parameter  larger  than  that  of  F  or  an 

crit 

air  content  larger  than  Tcr^ti  then  the  trajectories  corresponding 

to  a  Class  4  solution  would  be  found  outside  of  the  separatrix. 

Choosing  values  of  y  . ^  and  F_  near  zero  results  in  small 

'crit  C  .  _ 

crit 

scale  oscillations  about  the  vortex  point  and  is  classified  as  a 
Class  1  solution.  Critical  values  of  the  air  content  parameter 
and  forcing  function  parameter  evaluated  at  the  location  of  the 
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1 

1 

saddle  point  for  0  ■  0.30  are  Fr  ■  0.1834101  and  : 

crit  i 

Ycrit  ■  3.1022749.  Figure  22  is  a  plot  of  the  critical  parameters 

versus  the  normalized  bubble  radius.  In  this  study  we  are  J 

interested  in  values  of  u  >  1.  Therefore,  it  is  not  necessary  to  ^ 

begin  the  origin  of  u  at  zero.  It  is  important  to  note  that  the  ! 

magnitudes  calculated  for  the  critical  parameters  in  Fig.  22  and  the 

amplitude  and  frequency  of  oscillation  Indicates  in  Fig.  21  are  4 

» 

comparable  to  those  calculated  by  Parkin  [11].  This  apparent 
consistency  in  the  generalized  phase  plane  allows  us  to  consider  the 
more  specific  case  of  the  phase  plane  when  we  use  the  initial 
conditions  formulated  earlier.  Also,  if  we  let  u  +  »  then  we 

see  the  phase  plane  trajectories  asymptotically  approach  a  value  ! 

I 

equal  to  i 


u 


± 


h  Fc 

V3  n  *Q>2 


(3.39) 


The  asymptotes  are  represented  by  the  horizontal  dotted  lines  on 
the  phase  plane  plot  in  Fig.  21  and  apply  to  non-periodic 
trajectories  as  t  «. 

If  one  analyzes  the  phase  plane  when  the  initial  conditions  of 
Eq.  (3.11b)  are  used,  the  resulting  equation  is  similar  to 
Eq.  (3.30)  except  that  u(0)  ■  1.  The  phase  plane  equation  is 
written  as 


(3.40) 


CRIT 
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NORMALIZED  BUBBLE  RADIUS,  u 

Figure  22.  Critical  Parameters,  ycrit  and  FCrrit,  that  will  Produce  a  Bubble 
Growth  in  Accordance  with  the  Separatrix. 


which  satisfies  the  initial  conditions  exactly.  By  varying  the 

forcing  function  parameter  Fc  about  the  value  of  the  critical 

parameter  Fccrit*  one  obtains  the  phase  plane  plot  shown  in  Fig.  23. 

Notice  that  all  trajectories  begin  from  u  -  1,  u  -  0  and  follow  a 

trend  similar  to  the  trajectories  shown  in  the  general  phase  plane. 

As  noted  earlier  in  the  derivation  of  u(0),  the  value  of  u(0)  has 

_8 

a  magnitude  approximately  equal  to  1  x  10  fps.  The  fact  that  this 
initial  velocity  is  not  equal  to  zero  is  very  important  because  it 
allows  us  to  explain  the  existence  of  the  trajectories  when  there 
is  no  forcing.  If  one  looks  at  the  parabolic  coefficients  of  the 
forcing  function  represented  by  Eq.  (2.54),  it  is  easy  to  see  that 
the  zeroth  order  term  ■  0.0.  Therefore  if  the  bubble  has  zero 
wall  velocity  at  the  initial  point,  then  the  bubble  will  grow  only 
if  a  force  favorable  to  growth  is  applied.  Because  the  zeroth 
order  forcing  constant  is  zero  there  must  be  an  initial  wall 
velocity  enabling  the  bubble  to  traverse  along  any  one  trajectory 
as  shown  in  Fig.  23.  Otherwise  the  phase  plane  would  be 
represented  by  a  single  point  in  the  phase  plane  at  u  *  0,  u  *  1. 

Thus,  approximation  of  the  initial  velocity  as  zero  allows  one  to 
simplify  the  equations  without  sacrificing  the  accuracy  of  the 
calculation. 

In  summary,  we  have  looked  at  the  pha,?e  plane  plots,  potential 
energy  and  singularities  for  the  autonomous  zeroth  order  differen¬ 
tial  equation  Fq  >  0  and  F^  *  0.  The  singularities  which  exist  in  the 
system  for  >  0  explicitly  define  the  limits  between  periodic  and 
non-periodic  solutions.  The  character  of  these  singularities  was 
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Figure  23.  Phase-Plane  Trajectories  Originating  From  the 
Tnitial  Conditions. 
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verified  by  Liapunov's  Method  and  revealed  a  vortex  point  near  the  jj 

general  Initial  condition  u  -  u(0)  and  a  saddle  point  located  farther  & 

away  from  the  Initial  radius.  As  a  result,  the  phase  plane  plot  ! 

showed  trajectories  representing  periodic  solutions  circling  about 
the  vortex  point  while  the  saddle  point  acted  as  a  bridge  between  j 

the  periodic  and  non-periodic  trajectories.  Also,  the  location  j 

along  the  u-axis  of  the  vortex  and  saddle  points  was  exactly  the  J 

I 

t 

same  as  the  location  of  the  minimum  and  maximum  points  exhibited 
by  the  potential  energy  function.  Using  the  location  of  the 

saddle  point,  we  were  able  to  derive  from  the  potential  energy  i 

function  and  its  derivative  expressions  for  the  critical  parameters 
y  and  F_  which  are  used  in  Eq.  (3.40)  to  define  the 

c  C  .  _  | 

crit 

separatrix.  The  separatrix  which  is  now  a  function  of  Q  and  u  is 

used  to  show  the  separation  between  Class  1  and  Class  4  solutions.  J 

Thus,  the  trajectory  representing  the  separatrix  is  the  outermost 

limit  where  a  periodic  solution  exists.  All  solutions  within  the 

separatrix  are  periodic  motions  and  are  the  main  focus  of  this 

study.  And  all  trajectories  representing  periodic  or  non-periodic 

i 

i 

solutions  approach  an  asymptotic  value  indicating  a  constant  rate  1 

of  growth  or  collapse  of  a  bubble  as  its  radius  grows  infinitely.  ! 

i 

Having  established  limits  between  the  periodic  and  non¬ 
periodic  solutions,  one  is  now  ready  to  put  the  zeroth  order 
differential  equation  into  integral  form  in  order  to  determine 
an  approximate  solution  of  the  individual  bubble  oscillations. 
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3.2. 1.3  Polynomial  approximation  of  the  logarithmic  air  content 
parameter 

In  an  effort  to  solve  the  zeroth  order  differential  equation, 
one  oust  rewrite  Eq.  (3.40)  using  the  positive  root  only  as 

(u2  -  1 

3/2 
u 

(3.41) 

By  inverting  the  variables  in  Eq.  (3.41),  an  expression  for  the 
period  of  oscillation  can  be  written  in  the  form  of  a  definite 
integral  as 


) 


(1  +  0)' 


(u  -  0 


3/2. 
x  dx 


U’1 

V  (1  +  0) 


Jinx  - 


1 


(1  +  0) 


(x2  -  1)  +  4 


F 


(1  +  Q) 


C  f  3  .  \ 

2  (x  -  0 


(3.42) 


The  lower  limit  of  the  integral  is  the  initial  normalized  radius  while 
the  upper  limit  is  the  value  of  the  normalized  radius  u,  between 
u  ■  1  and  u  *  u8  where  ug  corresponds  to  the  maximum  value  of  u  when 
u  -  0  in  the  phase  plane  plot  in  Fig.  23.  As  it  stands,  Eq.  (3.42)  is 
not  soluable  by  any  standard  analytic  technique.  It  closely  resembles 
the  form  of  an  elliptic  integral  except  that  the  numerator  does  not 
have  a  whole  numbered  exponent  and  the  denominator  contains  a  natural 
log  term. 

In  order  to  shape  Eq .  (3.41)  into  a  form  which  may  be  solved  by 
some  known  analytic  technique,  one  must  first  multiply  the  numerator 
and  denominator  by  the  square  root  of  the  normalized  radius.  The 


integral  is  then 
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(3.43) 


which  still  contains  a  natural  logarithm  in  the  denominator. 
Generally,  one  would  numerically  integrate  Eq.  (3.43)  and  probably 
obtain  very  accurate  results.  Because  the  emphasis  of  this  study  is 
to  obtain  an  analytic  solution  to  the  Rayleigh-Plesset  equation, 
certain  modifications  of  Eq .  (3.43)  must  be  made  enabling  one  to  get 
an  analytic  solution  to  the  zero  order  equation.  If  one  considers 
the  function  u£nu,  which  is  contained  within  the  square  root  of  the 
denominator,  it  is  easily  seen  that  the  function  is  a  rather  smooth 
function  of  constantly  increasing  magnitude.  By  fitting  the  ufcnu 
function  with  a  cubic  polynomial  the  entire  denominator  can  be 
written  as  the  square  root  of  a  polynomial.  For  the  zero  order 
equation,  the  term  F^  -  0.  Then  the  polynomial  in  the  denominator 
would  be  a  cubic  . 

In  order  to  fit  the  u£nu  function  with  a  cubic  polynomial  four 
points  must  be  used  to  model  the  function  in  the  region  of  interest. 
Thus,  for  the  four  points,  uj,  U2»  u3*  u4»  one  writes 

3  2 

f  (uj)  -  UjinUj  ■  au^  +  bu^  +  cu^  +  d 

3  2 

f(u.)  -  u.£nu.  ■  au,  +  bu,  +  cu.  +  d 

L  L  (3.44) 

3  2 

f  (uj)  *  u^Jlnu^  *  au^  +  bu^  +  cu^  +  d 

3  2 

f(u4)  -  u4*nu4  -  au4  +  b«4  +  cu4  +  d  . 
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In  matrix  form,  Eq.  (3.44)  is  written  as 


jm 

" 

—  — 

—  - 

3 

2 

1 

u  ^  £nu  j 

U1 

U1 

U1 

a 

3 

U2 

2 

u2 
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1 

b 

u2lnu2 
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2 

1 

u^lnu  j 

U3 

U3 

U3 

c 

3 

u4 

2 

u4 

U4 

1 

d 

Vnu4 

— 

. 

»  . 

mm 

which  when  solved  with  the  IMSL  subroutine  called  LINV3F  produces 
the  coefficients 

a  -  -  0.050304 
b  -  0.578693 

c  -  -  0.001871 
d  -  -  0.526519  . 

The  numerical  values  of  uj,  U2»  uj,  u^  corresponded  to  the  values  of 
the  normalized  radius  starting  at  the  initial  condition  uj  ■  1  and 
ending  close  to  the  location  of  the  saddle  point  U4  -  3.60.  Points 
U2  and  U3  were  arbitrarily  chosen  to  be  U2  ■  1.80  and  U3  ■  2.60. 
Table  1  shows  the  %  error  between  the  function  uinu  and  the  cubic 
polynomial  f  ■  au^  +  bu^  +  cu  +  d.  It  is  evident  that  the  error 
falls  in  the  range  of  -1.0%  <  error  <  1.0%. 

To  see  the  effect  of  the  approximation  on  the  zeroth  order 
equation,  one  can  rewrite  Eq .  (3.40)  using  the  cubic  polynomial  to 
get 
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$ 
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TABLE  1 

Polynomial  Curve  Fit  of  u&nu 


u 

u&nu 

f  ■  au^  +  bu^  +  cu  +  d 

%  error 

1.0 

0.00000 

0.00000 

0.00000 

1.2 

0.21879 

0.21763 

-  0.53019 

1.4 

0.47106 

0.46707 

-  0.84794 

1.6 

0.75201 

0.74590 

-  0.81249 

1.8 

1.05802 

1.05171 

-  0.59681 

2.0 

1.38629 

1.38208 

-  0.30376 

2.4 

2.10112 

2.10686 

+  0.27318 

2.6 

2.48433 

2.49644 

+  0.48737 

2.8 

2.88293 

2.90092 

+  0.62408 

3.2 

3.72208 

3.74495 

+  0.61441 

3.6 

4.61136 

4.61963 

+  0.17927 
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Now,  Eq.  (3.46)  is  In  a  form  which  allows  one  to  calculate 
analytically  the  value  of  the  integral.  Solution  of  the  integral 
will  be  completed  in  the  next  section. 

3. 2. 1.4  An  Approximate  Zero  Order  Solution  using  Elliptic 
Functions 

Having  established  an  approximate  form  of  the  integral 
representing  the  period  of  oscillation  of  a  bubble,  an  analytic 
solution  in  the  form  of  elliptic  functions  must  be  formulated. 
Because  the  numerator  is  a  squared  term  and  the  denominator  is 
the  square  root  of  a  polynomial,  an  integral  from  Byrd  and 
Friedman  [131  can  be  applied.  The  solution  is  dependent  on 
the  roots  of  the  denominator  and  can  be  expressed  in  terms  of 
incomplete  elliptic  integrals,  the  first,  second  and  third  kinds 
and  a  product  of  Jacobian  elliptic  functions. 

The  first  thing  to  do  is  to  calculate  the  three  roots  of  the 
cubic  polynomial  in  the  denominator  of  the  integrand.  Because  the 
roots  can  vary  due  to  the  choice  of  y  and  Q,  we  will  choose 
Y  ■  1.4  and  Q  *  0.0  to  get  the  largest  possible  roots  available 
for  the  zero  order  solution. 

Using  Eq.  (3.46),  the  radicand  should  be  written  as 

3  2 

A[u  +  pu  +  qu  +  r ] 

where  p  ■  B/A,  q  -  C/A  and  r  *  D/A.  If  we  set 


which  denotes  the  phase-plane  location  of  the  vortex  point. 
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Then  the  coefficients  p,  q,  r  and  A  become 


2bu‘ 


P  - 


2au  -  1 
v 


q  ■ 


2cu  +  1 
v 


2au;  -  1 


2du2 


l 


(1  +  Q)' 


[2au^  -  1] 


Letting  the  three  real  roots  of  this  modified  cubic  be  u^,  and 
u^,  we  can  order  the  roots  as 


Ui  >  U  *  u2  >  u3 

with  u^  being  the  largest  root  and  being  the  smallest.  From 
Fig.  20  it  is  clear  to  see  that  all  the  trajectories  in  the  phase 
plane  originate  from  u  ■  1.  We  can  therefore  assume  that  the 
modified  cubic  always  has  a  root  at  unity.  Normally  the  cubic  can 
be  factored  using  Tartaglia's  Method,  but  because  we  know  one  of 
the  roots  is  u  ■  1  the  other  roots  can  be  found  by  dividing  the 
modified  cubic  by  u  -  1  and  solving  the  resulting  quadratic 


equation 
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Long  division  of  the  modified  cubic  by  u  -  1  gives 

x2  +  x(p  +  1)  +  (q  +  p  +  1) 
with  a  remainder  of 
2 

2u(a+b+c+d) 

2au2  -  1  ' 

V 

Then  with  the  numerical  values  of  the  coefficients  a,  b,  c  and  d 
from  the  curve  fit  of  uinu  (p.  85),  one  can  calculate  the  remainder 
to  be  equal  to  zero.  Therefore,  the  remaining  two  roots  are  found 
from  the  quadratic  equation.  Upon  factorizing  the  quadratic 
equation  and  applying  the  appropriate  coefficients,  the  remaining 
two  roots  are 

1  ± 

The  numerical  equivalents  of  these  two  roots  when  y  =  1.4  and 
Q  ■  0.0  are 

u  -  1.3662 

and 

u  -  -  0.9459  . 

Comparison  of  all  three  roots  with  the  range  of  integration  as 
prescribed  by  the  phase  plane  of  Fig.  20  justifies  the  ordering 
of  the  roots  as 


1.0568u2  -  1 
_ v _ 

2(1  +  0. 100608u2) 


4.212144u2(l  +  0.10O608u2) 
v'  v' 


( 1.0568u^  -  1) 


Uj  =  1.3662  ,  U£  =  1.00  ,  u^  =  -  0.9459 


« 

< 
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Consequently  Eq.  (3.46)  can  be  rewritten  using  u^,  u^  and  u^  as 


T  m  ->  /  (1  ±  Q)3  Ur  _ 

V  (2au^  -  1)  u_!  (xl  " 


2. 

x  du 


x)(x  -  x 2 ) (x  -  x3) 


(3.48) 


Since  all  trajectories  shown  in  Fig.  20  originate  at  u  ■  1,  the 
analysis  conducted  above  for  the  trajectories  that  lie  to  the 
right  of  u  »  1  should  also  work  for  the  trajectories  that  lie  to 
the  left  of  u  ■  1. 

It  has  been  previously  stated  (p.  75)  that  the  trajectories 
to  the  left  of  u  -  1  represent  small  scale  oscillations  of  a 
flaccid  air  bubble.  Indeed,  the  main  focus  of  this  work  is  to 
study  the  bubble  growth  as  it  occurs  in  the  vaporous  growth 
region.  However,  in  an  effort  to  maintain  completeness  in  the 
solution  of  the  zero-order  equation  one  should  be  able  to  use  a 
similar  factorization  as  shown  above  to  define  the  ordering  of 
the  roots  for  the  trajectories  to  the  left  of  u  ■  1.  These 
roots  could  then  be  used  in  the  solution  of  Eq.  (3.48)  to 
describe  the  time  history  of  these  small  scale  flaccid  air 
bubble  oscillations. 

Considering  the  same  factorization  as  before,  the  roots  are 
written  as 


u:  -  1.0 
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2,3 


2au2  +  2bu2  -  1 
v _ v 

2(1  -  2an2) 


1  i  \/  1  + 


8u2(a  +  b  +  c)  -  16u^(a2  +  ab  +  ac) 
v  v 


( 2au2  +  2bu2  -  l)2 


3  3 

Using  a  polynomial  curve  fit  for  the  function  u£nu  *  au  +  bu  +  cu  +  d 


where 

a  -  -  0.205484 
b  -  1.111631 

c  -  -  0.606996 
d  -  -  0.299157 

over  the  range  u  *  0.8  to  u  ■  1.0,  one  can  write 
1 . 8 1 2  3u  2  -  1 

V 

q  *  o 

,-5  2(1  +  0.41095u^) 


2.  39326u2(  1  +  0.41095u2) 
Vv  vy 

(•  8123u2  -  l)2 


Choosing  y  ■  1*4  and  0  *  0.30  then  one  finds  the  parameter  u^  equal  to 

u  *  0.910166  and  the  roots  are 
v 


Uj  =  1.00 


u  2  =  0.823184 


u 


3 


-  0.449188 


Therefore,  the  ordering  of  the  roots  as 


04 


Uj  >  u  > 

is  equivalent  to  the  ordering  of  the  roots  for  the  traiectories 
that  lie  to  the  right  of  u  *  1.  Comparison  of  the  roots  u^, 
and  Uj  against  the  roots  as  calculated  by  numerical  methods 
resulted  in  errors  of  approximately  0.60%.  Thus,  the  factorization 
method  is  a  reliable  means  by  which  the  roots  of  the  cubic 
polynomial  in  the  denominator  of  Eq .  (3.46)  can  be  calculated. 

It  is  interesting  to  note  that  when  values  of  y  and  Q  are 
chosen  such  that  the  parameter  uv  increases  in  magnitude  from  less 
than  one  to  greater  than  one,  the  roots  u^,  U£  and  u^  migrate  from 
left  to  right  along  the  u-axis.  When  the  parameter  uv  equals  unity 
there  is  a  double  root  at  u  ■  1  which  sets  the  upper  and  lower 
limit  of  Eq.  (3.48)  equal  to  unity.  Previous  analysis  reminds  us 
that  for  a  double  root  at  unity  a  bubble  radius  will  neither 
increase  nor  decrease  in  size  as  time  increases.  Thus,  it  will  be 
at  this  point  that  the  analysis  will  continue,  considering  only 
values  of  the  parameter  uv  greater  than  or  equal  to  unity,  even 
though  an  equivalent  analysis  could  be  performed  for  values  of  uv 
less  than  one.  Refering  to  Byrd  and  Friedman,  Eq.  (3.48)  can  be 
written  as 


a 
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where 

V  -  -75-  [E(*,k)  -  k2SnuCdu] 

1  k  1 

(3. 49) 

V  -  “TT  [2(2  -  k2)E(^,k)  -  k'2F(*,k) 

1  3k  * 


-  k2SnuCdu(k  2nd2u  +  4  -  2k2)] 

The  parameters  in  Eq.  (3.49)  are  functions  of  the  roots,  uj, 
U2  and  U3  and  are  defined  as  follows: 


Modulus 


2 


Note:  u  is  the  upper  limit  of  integration. 


¥VWV*Virv*’U  ru  »V«'V  wmkw^j  m\j  w"w  ru  p 


1  M* 


\ rw v r.TWW  W  H/l  *n  Kn  *  *cn ¥ n in  *  * 


T  -  F($,k)  [gu?  -  gu2(l  -  -~] 

J  L  2  3k  z 


+  E(*,k)[2gu  u  (l  -  *  gu^(l  -  (-2-(2  hr-1)] 

1  J  u2  k  z  z  u2  3k  4 


-  SnuCdu[2gu  u  (l  -  — -H^y)  +  gu^l  -  ~)  (“T^X4  ~  2k2)] 
Li  u2  k  z  L  u2  3k  4 


9  9  ^  9  V  ^ 

SnuCdund  u[gu,(l  -  —■)  ( — pr)] 

2  3k  L 


(3.50) 


Now,  if  one  considers  the  terms  inside  the  square  parentheses 
in  Eq.  (3.50)  for  a  particular  value  of  Q  and  y»  the  quantities 
inside  the  square  parentheses  are  considered  constant.  All 
combinations  of  the  elliptic  functions  in  front  of  the  square 
parentheses  are  functions  of  $  which  are  functions  of  the  variable 
u.  Therefore,  x  ■  f(u)  which  is  the  inverse  of  what  we  want;  namely, 
u  ■  f(x).  Also,  since  all  the  elliptic  functions  in  Eq.  (3.50)  do 
not  have  well  defined  inverses,  an  alternate  strategy  must  be  used 
to  invert  Eq.  (3.50)  in  order  to  write  an  expression  describing  the 
growth  of  the  normalized  bubble  radius  versus  the  bubble  time. 

Since  all  terms  in  the  square  parentheses  in  Eq.  (3.50) 
represent  constants,  one  can  rewrite  Eq.  (3.50)  as 


X  *  F($,k)V^  +  E($,k)V2  ~  SnuCduV^  -  SnuCdund  uV^ 


(3.51) 
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where  V^,  and  are  constants  for  a  particular  value  of  0 


and  y.  Now  since  the  sum  of  the  elliptic  functions  in  Eq.  (3.51) 
does  not  have  a  known  inverse,  an  approximation  of  each  function 
must  be  derived  that  allows  one  to  invert  the  equation  to  get 
u  «f(r).  Because  the  functions  F($,k)  and  E(<fi,k)  are  defined  as 
definite  integrals,  a  numerical  routine  must  be  used  to  evaluate 
the  functions  across  the  range  of  the  roots  uj  and  u£.  The 
functions  F(<ji,k)  and  E($,k)  are  considered  incomplete  elliptic 
integrals  when  $  <  ir/2.  If  $  ■  ir/2  then  the  functions  are 
considered  complete  and  are  tabulated  in  Byrd  and  Friedman's 
Table  of  Elliptic  Integrals  (13].  Unfortunately,  the  tables  are 
not  comprehensive  enough  to  Include  every  possible  value  of 
F(<j>,k)  and  E($,k)  for  different  values  of  <f>  so  that  the  numerical 
routine  utilizing  Simpson's  1/3  Rule  was  used. 

In  Fig.  25,  the  function  F($,k)  was  plotted  as  a  function  of 
the  modular  angle  $.  Because  the  values  of  Q  and  y  vary  in 
accordance  with  the  different  phase  plane  trajectories  located  to 
the  right  of  the  initial  condition  as  shown  in  Fig.  20,  extreme 
values  of  k  were  computed  to  represent  those  trajectories.  Using 
these  extreme  values  of  k,  curves  representing  the  maximum  and 
minimum  range  that  the  function  F($,k)  would  fall  are  plotted  as 
solid  lines  in  Fig.  25.  All  other  values  that  F(<j>,k)  might  attain 
would  lie  in  between  those  two  solid  lines.  Recall  the  modulus  k 
is  defined  by 


Vui  -u 


(3.52) 


K 


'  ^  ‘»V  *.*  * 


where  u^,  and  are  the  roots  of  the  cubic  polynomial  in  the 
denomirator  of  Eq.  (3.46).  When  Q  ■  0.183216  and  y  ■  1.4,  the 
minimum  potential  energy  falls  directly  on  the  initial  condition. 

The  cubic  polynomial  in  the  denominator  of  Eq .  (3.46)  has  double 
root  at  uj  ■  U2  ■  1.0.  Consequently,  k  ■  0  is  the  minimum  value  of 
k  used  to  plot  the  extreme  upper  value  of  the  elliptic  function. 
Likewise,  when  Q  -  0.0  and  y  *  1.4,  the  cubic  polynomial  has  roots 
such  that  u^  >  U2  >  U3.  The  result  is  a  maximum  value  of 
k  -  0.446046  which  is  used  to  compute  the  extreme  lower  value  of  the 
elliptic  functions.  This  same  procedure  can  be  used  for  the  other 
elliptic  functions  that  come  out  of  Eq.  (3.51).  Notice  that  for 
F($,k)  and  E(<j>,k),  the  maximum  curve  corresponds  to  k  *  0.  The 
opposite  is  true  when  looking  at  the  curves  for  SnuCdu  and 
SnuCdund^u  in  Figs.  27  and  28. 

Having  determined  a  maximum  and  minimum  range  for  the  elliptic 
functions,  it  was  of  great  benefit  to  derive  a  separate  functional 
representation  for  each  elliptic  function  shown  in  Figs.  25  through 
28.  To  do  this,  a  median  curve  lying  in  the  middle  of  the  two  solid 
lines  of  each  figure  was  plotted.  The  median  curve  is  shown  as  a 
dotted  line  and  was  calculated  using  k  ■  0.315759  as  the  modulus. 
Each  median  line  is  used  to  represent  the  entire  range  of  each 
elliptic  function  since  there  could  be  an  infinite  number  of 
representative  curves  for  all  combinations  of  Q  and  y  chosen. 
Finally,  since  each  median  line  is  roughly  parabolic  in  shape,  a 
parabolic  approximation  of  each  dotted  line  in  Figs.  25  through  28 
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was  derived.  The  coefficients  were  derived  by  fitting  each  dotted 
line  at  three  representative  points  starting  with  the  zero  point. 
Using  the  following  parabolic  equations 

F(*,k)  -  A^2  +  +  CL 


E( 4» ,k)  -  A2$  +  B24>  +  C2 


SnuCdu  ■  A^<^  +  B^ifr  + 


2  2 

SnuCdund  u  *  A^$  +  B^4»  + 


the  coefficients  were  computed  as 


At  -  0.021844 


B  -  0.992116 


A  -  -  0.02073 


B2  -  1.007155 


A  -  -  0.83226 


B3  -  1.307305 


A  -  -  0.88065 
4 


B.  *  1.38332 

4 


C1  ■  C2  ’  C3  ■  C4  ■  °-° 


(3.53) 
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Figures  29  through  32  are  plots  of  the  resultant  curve  fits 
(dotted  line)  of  the  median  representative  of  the  various  elliptic 
functions  (solid  lines).  Now,  one  can  rewrite  Eq.  (3.51)  as 

T  -  V^Aj*2  +  +  Cj)  +  V2(A2*2  +  B2*  +  C2) 

+  V3(A3^2  +  B3^  +  C3)  +  V4(A4*2  +  Ba«^  +  C4)  (3.54) 

which  has  no  term  larger  than  second  order  in  4>. 

Grouping  like  powers  of  $,  we  write  Eq.  (3.54)  in  the  general 
form  of  a  parabola  as 

t  *  a$2  +  b<ji  +  c  (3.55) 


where 


a  -  Vl  +  V2A2  -  V3A3  -  V4A4 
b  -  VjBj  +  V2B2  -  V3B3  -  V4B4 


c  «  VjCl  +  v2c2  -  v3c3  -  v4c4  . 


It  should  be  noted  that  Eq .  (3.55)  is  now  written  as  a  function  of 
p.  But,  because  the  definition  of  the  modular  angle  $  is  written 
as  a  function  of  u,  we  still  have  t  ”  f(u).  It  is  desired  to  get 
the  inverse  f  -  g(i).  Therefore,  Eq .  (3.55)  can  be  inverted  by 
using  the  quadratic  formula.  Subtracting  t  and  applying  the 
quadratic  formula,  one  gets 


* 


-  b  ±  V  b2  -  4a(^  -  t) 

2a 


(3.56) 
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no 

Since  -  0  and  $  must  equal  zero  at  t  ■  0,  one 

must  choose  the  positive  branch  of  the  square  root  in  Eq.  (3.56)  to 
represent  the  inverse.  Having  established  $  ■  f(x),  one  can  now  use 
the  fact  that 


sin  $ 


(3.57) 


to  derive  an  expression  in  the  form  u  -  f(x).  Solving  for  u  from 
Eq.  (3.57),  one  gets 


u 


K-uJ 

L  9  9  ]  +  U3 


2  2 
k  sin  *  -  1 


(3.58) 


Substituting  expressions  for  k(uj,u2,u3)  and  one  can  write 
u  ■  f(x)  as 


u3  ”  u2 


£4^)  »m2 


+  u. 


2a 


(3.59) 


where  uj,  U2  and  U3  are  functions  of  Q,  Fq  and  y  as  defined  by 
Eq.  (3.25a).  The  resultant  plot  of  the  normalized  bubble  radius  as 
a  function  of  the  bubble  time  is  shown  in  Fig.  33.  Each  curve 
beginning  at  x  *  0  corresponds  to  a  particular  trajectory  in  the 
phase  plane.  Each  trajectory  is  characterized  by  the  author's 
choice  of  Q  and  y  in  such  a  manner  that  the  trajectory  is  always 
located  at  or  to  the  right  of  the  initial  condition.  Since  each 
trajectory  is  centered  about  a  vortex  point  whose  location  along 
the  u-axis  is  defined  by 

u  =  -  , 

v  (1  +  Q) 


DIMENSIONLESS  BUBBlf  TIME,  T 

Figure  33.  A  Plot  of  Numerical  Data  Showing  the  Zero  Order 
Solution,  u  ”  u(T) . 
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each  curve  in  Fig.  33  has  a  particular  value  of  uv  labeled  at  the 
half  period.  If  the  calculation  had  been  carried  across  the  upper 
and  lower  half  of  the  phase  plane  trajectory,  then  all  the  curves  in 
Fig.  33  would  return  to  the  x-axis  completing  a  whole  period  of 
oscillation. 

It  is  interesting  to  note  that  these  oscillations  are  of  a  very 
high  frequency  and  as  a  bubble  passes  through  a  region  favorable  to 
vaporous  growth  or  collapse,  these  free  oscillations  occur  a  number 
of  times  before  the  bubble  reaches  the  separation  point.  From 
Fig.  11,  one  can  determine  that  at  K  *  0.60  the  dimensionless  time 
from  the  initial  point  up  to  the  separation  point  equals 
approximately  760  units  of  bubble  time,  i.e.,  x  -  760.  If  the 
complete  period  of  oscillation  is  Tperi0(j  x  6.2  units  of  bubble 
time,  one  could  calculate  the  number  of  individual  free  oscillations 
to  be  NCyCies  ■  122.  In  laboratory  time,  this  means  that  a 
particular  bubble  will  oscillate  *  122  times  over  a  time  period  of 
0.144  seconds. 

Having  established  some  knowledge  about  the  zeroth  order 
solution  in  terms  of  elliptic  functions,  one  now  desires  a  more 
usable  form  of  the  solution,  preferably  in  terms  of  trigonometric 
functions,  that  will  closely  approximate  Fig.  33.  Looking  at  the 
curves  in  Fig.  33,  one  could  say  they  are  roughly  of  the  shape 
u  -  1  -  cos  x.  Using  a  more  general  form  of  the  curve  1  -  Cos  x, 
one  can  write 

u  ■  A  +  B  cos  (it  x/ xf  )  +  C  cos(2ir  x/xf)  (3.60) 
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where  is  the  dimensionless  time  corresponding  to  the  time  it 
takes  to  reach  the  end  of  a  half  period  and  A,  B  and  C  are 
coefficients  to  be  determined.  Using  the  following  three  conditions 

@  t  ■  0  ,  u  -  1  , 


(?  T  ■  t  ,  u  ■  u  , 
m  m 


and 


@  t  -  xf  ,  u  -  uf  , 


one  can  solve  simultaneously  for  A,  B  and  C  to  get 


A  -  |  (!  +  uf )  - 


u  - 
m 


{  ([  +  Uf)  -  }  (!  “  Uf)C08("  Tn/Tf  ) 

cos  (2ir  r  /re)  -  1 
'  mi 


B  -  J  (1  “  uf)  (3.61) 

C  U°>  "I  (1  +  Uf|  ~  I  ^  ~  uf)C0S^  Tm/Tf) 

COS  (2 7T  T  lie  ]  -  1 

'  mi' 

The  parameters  um  and  Uf  correspond  to  the  values  of  u  alone  a  line 
representing  midpoint  values  of  the  zeroth  order  solution  and  the 
values  of  u  along  a  line  representing  values  at  the  end  of  a  half 
period.  Figure  34  shows  the  lines  along  which  the  values  of  ura  and 
Uf  were  taken.  Upon  plotting  these  values  of  um  and  Uf  against  the 
parameter  uv,  the  following  curve  fits  were  derived: 


u  =  -  4.0477  u2  +  9.2935  u  -  4.2458 
ra  v  v 


u  =  1.99873  u  -  0.99873  . 

f  v 


(3.62) 

(3.63) 


Stf 
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The  times  that  correspond  to  the  values  um  and  uf  are  defined  as 
Tjg  and  Tf  and  are  read  off  the  same  lines  that  um  and  Uf  were  in 
Fig.  34.  Plotting  Tm  and  Tf  against  the  parameter  uv  resulted  in  a 
parabolic  curve  fit  for  Tm  as 

t  ■  1.25  +  £  sin  0  (3.64) 

m 

where 

£  -  -  4.3239(uy  -  l)2  +  1.2804(uv  -  l) 
and  a  linear  fit  for  Tf  as 

T,  -  2.1587  +  5.6764  u  sin  6  .  (3.65) 

f  v 

The  parameter  0  defined  the  angle  that  the  line  along  which  uf  and 
Tf  deviated  from  the  vertical.  Measurement  of  0  gives  0  ■  1.208 
radians.  Figures  35  and  36  show  the  approximate  curve  fits  of  the 
parameters  um  versus  uv  and  £  versus  (uv  -  l).  The  curve  fits  show 
relatively  good  accuracy  and  thus  can  be  used  with  confidence  in 
determining  an  approximation  to  the  zeroth  order  solution. 

Therefore,  plotting  Eq.  (3.60)  in  combination  with  Eqs.  (3.61), 
(3.62),  (3.63),  (3.64)  and  (3.65)  against  the  dimensionless  bubble 
time  t  gives  a  maximum  error  of  -  4.5%  of  the  zero  order  solution 
as  shown  in  Fig.  37.  In  order  to  reduce  the  error  even  more, 

Eq.  (3.60)  was  rewritten  as 

u  ■  (l  -  e)[A  +  B  cos(ir  t/t^)  +  C  cos(2w  t/t^)]  (3.66) 

where  e  is  some  functional  approximation  of  the  relative  error  used 
to  reduce  the  error  between  the  solid  and  dotted  lines  shown  in 


Fig.  37.  Writing  e  as 
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e  ■  A  +  B  cos(ir  t/ )  +  C  cos(2ir  t/t^)  (3.67) 

and  solving  for  the  coefficients  a  new  approximation  of  the  zeroth 
order  solution  can  be  written.  Using  the  conditions  on  e  that 

e  ■  0  @  t  ■  0 

e  -  e  @  t  ■  t 
m  ra 

e  -  0  @  t  -  rf 

the  coefficients  A,  B  and  C  are  written  as 

e 

.  _  _ m _ 

1  -  cos  (2ir  T  /  T,  1 

B  -  0  (3.68) 

e 

c  _ _ m 

1  -  cos  f 2ir  x  / Te  J 
v  mi'' 

where  em  is  the  maximum  error  indicated  for  each  value  of  uv  used. 
Plotting  the  values  of  em  versus  uv,  one  can  find  an  approximate 
curve  fit  for  em  which  is  written  as 

e  -  0.09774  u2  +  0.1585  u  -  0.1136  .  (3.69) 

m  v  v 

The  parameters  Tm  and  tf  were  defined  previously  by  Eqs.  (3.64)  and 
(3.65).  A  plot  of  Eq.  (3.66)  with  (3.67),  (3.68)  and  (3.69)  gives 
an  approximation  of  the  zeroth  order  solution  with  a  maximum  error 
of  -  1.5%  for  all  possible  values  of  uv.  Figure  38  shows  the 
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approximate  and  correct  curves  of  the  zeroth  order  solution  plotted 
as  functions  of  the  bubble  time  t.  Thus,  a  fairly  accurate  closed 
form  solution  exists  which  can  now  be  used  to  analyze  the  first 


order  equation. 
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Having  found  a  partial  multiple  scales  solution  of  the  pulse 

forced  Raylelgh-Plesset  equation  of  cavitation  bubble  dynamics,  we 

must  discuss  and  summarize  the  features  of  this  work.  First,  before 

any  analysis  was  performed,  certain  assumptions  were  made.  It  was 

assumed  that  one  can  model  the  dynamics  of  bubble-ring  cavitation  as 

though  each  bubble  were  a  small  isolated  spherical  bubble 

distributed  throughout  the  flow.  For  the  microbubbles  of  radius, 

Rq,  distributed  in  the  flow,  all  interactions  with  adjacent  bubbles 

and  walls  were  neglected.  This  assumption  was  made  strictly  in 

the  interest  of  simplifying  the  problem.  Also,  it  is  assumed 

cavitation  is  initiated  from  these  microbubbles  of  radius,  R  ,  in 

o 

the  minimum  pressure  region  which  is  on  the  boundary  of  the 
submerged  body.  This  fact  was  reflected  in  the  formulation  of  the 
dynamical  equation  by  using  the  pressure  distribution.  Recent 
measurements  made  by  Holl  and  Carroll  [1]  revealed  a  very  low 
pressure  area  along  the  surface  of  the  body  which,  under  special 
conditions,  created  an  environment  favorable  for  bubble 
growth.  Using  the  specific  pressure  coefficient  data  collected  by 
Holl  and  Carroll  [1],  a  forcing  function  was  derived  as  a  function 
of  the  dimensionless  bubble  time,  t.  The  forcing  function,  F(r), 
produced  a  pulse-like  force  which,  when  incorporated  into  the 
problem,  caused  the  dynamical  equation  to  be  non-autonomous .  Past 


studies  by  Parkin  [11]  were  made  using  a  piecewise  autonomous  step 
function  to  approximate  the  forcing  function.  This  alleviated  some 
of  the  problems  that  would  have  arisen  if  the  forcing  function  were 
modeled  as  a  parabola,  which  more  closely  represents  the  true  nature 
of  the  pulsed  force.  In  an  effort  to  capture  the  parabolic  nature 
of  the  forcing  function,  a  parabolic  fit  was  derived  which  satisfied 
the  physical  requirements  that  the  forcing  function  exhibited  in  the 
region  where  vaporous  growth  occurred. 

To  this  point,  the  development  of  the  problem  is  straight¬ 
forward.  A  differential  equation  and  initial  conditions  exist 
that  describes  the  dynamics  of  growth  along  with  a  forcing 
function  representing  the  pulsed -force  that  acts  on  the 
microbubbles.  The  problem  seems  like  an  ideal  candidate  for  a 
computer  solution  where  detailed  histories  of  the  bubble  growth 
could  be  plotted  using  an  appropriate  numerical  routine.  The 
problem  with  a  numerical  solution  is  the  loss  of  the  parametric 
analysis  of  the  critical  parameters  that  govern  cavitation 
inception  in  this  type  of  flow.  Also,  the  differential  operator 
of  the  Rayleigh-Plesset  equation  does  not  directly  lend  itself 
to  a  nice  solution  in  terms  of  common  well-behaved  functions. 
Therefore,  parametric  analysis  of  the  various  solutions  coming 
from  the  governing  dynamical  equation  is  nearly  impossible. 

Figure  39  is  a  prime  example  of  how  the  solution  varies  with 
respect  to  the  various  flow  parameters.  The  curves  in  Fig.  39 
represent  isothermal  cavitation  bubble  growth  in  response  to 


single  parabolic-pulse  forcing  functions.  The  transient  response 
is  directly  affected  by  the  relative  pulse  height  and  relative 
pulse  duration.  By  understanding  the  parametric  relationships 
that  define  the  growth,  one  can  better  analyze  the  specific  flow 
characteristics  without  relying  on  the  computer  to  produce 
detailed  histories  of  growth. 

On  the  other  hand,  certain  information  can  be  derived  from 
the  numerical  solutions  as  shown  in  Fig.  39.  The  curves  show 
essentially  no  growth  for  certain  values  of  the  parameter  8  while 
higher  values  of  8  cause  explosive  growth.  One  can  predict  the 
existence  of  two  different  classes  of  bubble  growth  where  the 
interface  between  the  two  classes  of  growth  is  a  particular  cutoff 
value  of  8.  Analysis  of  the  governing  equations  should  then 
reveal  a  parametric  relationship  that  allows  one  to  predict  that 
cutoff  value  of  8*  It  should  be  noted  that  the  curves  in  Fig.  39 
are  unpublished  results  obtained  by  Parkin  under  the  direction  of 
M.  S.  Plesset.  Also,  the  parameters  shown  in  Fig.  39  pertain  to 
Parkin's  work. 

Using  a  multiple  scales  analysis,  a  series  of  differential 
equations  and  initial  conditions  representing  the  zeroth,  first 
and  second  orders  of  e  are  derived.  From  the  zeroth  order 
equation,  the  global  representation  of  the  phase-plane  trajectories 
and  corresponding  vortex  and  saddle  points  are  found.  If  the 
autonomous  form  of  the  trajectory  equation  is  plotted,  then  the 
saddle  point  disappears  and  the  trajectories  become  closed-loop 


MW.V.W  Rl  M.WIWWWVRW.  V\  I^^WVWWnniOTUWWWYVruroTrv.rii  irjxu 


127 

paths  centered  around  a  vortex  point.  The  location  of  this  vortex 
point  can  be  located  rather  precisely  using  Liapunov's  method  or 
the  corresponding  potential  energy  curve. 

An  important  result  of  the  phase-plane  analysis  for  the 
autonomous  zeroth  order  equation  is  that  a  critical  condition 
was  found  which  represents  the  point  of  separation  between  two 
physically  distinct  types  of  motion.  The  trajectories  that  lie 
to  the  left  of  this  critical  point  are  representative  of  small 
scale  flaccid  air  bubble  oscillations  while  those  to  the  right 
have  a  somewhat  larger  amplitude  of  oscillation  due  to  vaporous 
growth.  The  critical  parameter  that  defines  the  various 
trajectories  is 


^TL 
1  +  Q 


(3.31) 


In  particular,  the  parameter,  uv ,  defines  the  location  of  the 
vortex  point  for  each  trajectory.  For  uy  greater  than  one,  the 
vortex  point  and  corresponding  trajectory  lie  to  the  right  of  the 
initial  condition.  For  uv  less  than  one,  the  vortex  point  and 
trajectory  lie  to  the  left  of  the  initial  condition.  Because  this 
work  concerns  itself  primarily  with  vaporous  growth,  we  have 
considered  only  values  of  greater  than  one. 

Having  defined  the  phase-plane  trajectories  for  the  zeroth 
order  autonomous  form  of  the  Rayleigh-Plesset  equation,  an 
approximate  solution  to  the  zeroth  order  equation  was  found.  From 
the  trajaectory  equations  one  can  write  the  bubble  time  equation 
for  t  in  quadrature  form  as 
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r 


2. 

x  dx 


u-1  / _ l 

V  (1  + 


xlnx  - 


1 


o) 


(1  +  0)' 


f  3  \ 

-  x) 


(4.1) 


Equation  (4.1)  is  generally  solved  numerically,  but  since  a 
parametric  solution  is  desired,  an  approximate  solution  is  obtained 
by  fitting  the  uJtnu  function  with  a  cubic  polynomial.  The  resulting 
integral  is  solved  using  a  Table  of  Elliptic  Integrals.  The 
solution,  which  is  expressed  in  terms  of  incomplete  elliptic 
integrals  of  the  first,  second  and  third  kinds  and  a  product  of 
Jacobian  elliptic  functions,  is  dependent  on  the  roots  of  the 
denominator  of  Eq.  (4.1).  The  roots  of  the  denominator  are  found 
rather  easily  using  the  fact  that  one  of  the  roots  is  equal  to 
unity.  This  stems  from  the  phase-plane  trajectories  which  all 
originate  from  the  point  u  ■  1.  Thus,  the  modified  cubic  in  the 
denominator  can  be  reduced  to  a  quadratic  form  which  is 
subsequently  factored  using  the  quadratic-equation.  The  roots 
are  then  ordered  from  the  largest  to  the  smallest  which  allows 
one  to  calculate  the  appropriate  solution  to  the  integral. 

Since  the  solution  is  not  readily  invertible,  an  approximate 
inverse  solution  of  the  form 


u-(l-e)[A  +  B  cos  ut  +  C  cos  2wt ) 


(4.2) 


where  A,  B,  C  and  w  are  functions  of  uy.  The  error  parameter,  e, 
simply  reduces  the  error  of  the  approximate  trigonometric  inverse. 
Use  of  Eq.  (4.2)  allows  one  to  look  at  one-half  of  a  period  for 
a  particular  value  of  uv<  Thus,  a  parametric  study  of  the 


msmmsfmm  wwwwwwuuiuwu 


WWW.  IMAflfflW  1/iVW  w^wwwi/*  vn  \rw  v*>r*  *  r-*  ™  nxwwrw  «ru  r-v  »-w  ■  »▼ 


129 

approximate  solution  u  as  a  function  of  the  parameter  u^.  Therefore 
an  approximate  solution  to  the  zeroth-order  autonomous  Rayleigh- 
Plesset  equation  exists  that  describes  the  growth  of  the  bubble 
radius  as  a  function  of  bubble  time,  t. 

It  is  important  to  note  that  this  work  concerns  itself  solely 
with  cavitation  inception.  At  no  time  does  this  formulation  allow 
for  recirculation  of  the  nuclei  within  the  laminar  separation 
bubble.  If  the  theory  were  concerned  with  desinence,  then  the 
dynamical  equation  and  solution  would  be  very  different. 

Finally,  after  considering  the  approximations  that  have  been 
made  to  derive  a  solution  to  the  zeroth  order  Rayleigh-Plesset 
equation,  a  comparison  of  the  results  is  made.  One  can  see  there 
is  good  agreement  between  the  amplitudes  of  oscillation  shown  in 
Fig.  39  with  the  amplitudes  shown  in  Fig.  33.  This  apparent 
consistency  between  the  older  results  of  Fig.  39  and  the  results 
achieved  by  use  of  the  multiple  scales  analysis  is  encouraging. 
Though  the  results  of  Fig.  33  do  not  show  the  explosive  bubble 
growth  that  the  older  results  exhibit,  this  matter  will  be  dealt 


with  in  the  first-order  solution 
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APPENDIX  A.  FORTRAN  CODE  DISCUSSION 


A1 .  The  Basic  Algorithm 

In  the  basic  formulation  of  the  problem,  it  was  required  that 
we  could  perform  a  parabolic  curve  fit  as  well  as  numerically 
integrate  the  experimental  Cp  data  obtained  by  Holl  and  Carroll 
[1].  The  form  of  the  integral  arose  when  an  expression  for  the 
dimensionless  bubble  time  was  derived  using  the  kinematics  of  the 
flow.  During  the  development  of  the  Fortran  code,  care  was  taken 
to  satisfy  all  the  criteria  laid  down  in  Chapter  2.  The  code  was 
designed  to  handle  an  odd  or  even  number  of  data  points  and  to 
output  values  of  the  computed  integral,  the  coefficients  of  the 
parabolic  curve  fit,  the  integrand  and  the  experimental  data. 

The  first  step  of  the  code  read  in  the  data  in  a  form 
x^,  f(x^).  These  data  were  then  used  to  compute  the  coefficients 
A,  B  and  C  of  the  curve  fit.  The  coefficients  were  then  used  to 
calculate  the  integral  contained  within  the  bubble  time  equation, 
Eq.  (2.8). 

Calculation  of  the  coefficients  A,  B  and  C  was  performed  by 
expanding  a  third  order  Lagrange  Polynomial,  Pj(s),  where 


tr  (s  -  s  ) 
i-0 


(A. 1.1) 


w  (s  -  s  ) 
1-0  J 


i*0 
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The  polynomial  was  then  multiplied  by  the  functional  value  f(s.j) 
for  each  data  point  and  summed  n  times,  where  n  corresponds  to  the 
degree  of  the  polynomial  approximation  desired.  Thus,  the 
polynomial  has  the  form 
n 

p  (s)  -  Z  f (s  )P  (s)  (A. 1.2) 

n  ,  —  J  J 

j-o 

where  p  (s)  is  the  desired  polynomial  of  degree  n  which  exactly 
n 

fits  n  data  points.  Equation  (A. 1.2)  represents  the  form  of  the 
polynomial  used  to  fit  the  experimetal  data  taken  by  Holl  and 
Carroll  [1].  The  formulae  for  these  calculations  are  derived  by 
expanding  Eq.  (A.  1.2)  and  using  n  ■  3  to  perform  a  parabolic 
approximation  across  n  -  3  data  points.  A  parabolic  fit  was  derived 
for  each  set  of  the  three  consecutive  data  points  beginning  with  the 
first  data  point.  The  coefficients  A,  B  and  C  were  calculated  for 
consecutive  triads  of  data  points  until  all  points  were  used.  The 
result  is  a  set  of  coefficients  A,  B  and  C  that  represent  the 
coefficients  of  a  piecewise  parabolic  approximation  of  consecutive 
triads  of  data.  This  type  of  curve  fit  proved  advantageous  for 
fitting  a  curve  over  the  irregular  portions  of  the  pressure 
distribution  while  use  of  a  single  parabola  resulted  in  a  poor  fit 
due  to  its  inability  to  capture  the  irregularities  near  the 
separation  point  and  the  pressure  recovery  region. 

Since  a  smooth  representation  of  the  pressure  distribution 
along  the  arc  of  the  body  was  desired,  several  points  were 
interpolated  between  the  experimental  points  using  the  appropriate 


set  of  parabolic  coefficients  derived  from  the  curve  fit.  The 


resulting  curve  fit  pressure  distribution  along  with  the  experi¬ 
mental  data  taken  by  Holl  and  Carroll  [1]  are  shown  in  Fig.  2.  By 
using  the  interpolated  pressure  distribution,  a  better  represent¬ 
ation  of  the  initial  conditions  was  made. 

Having  calculated  a  piecewise  parabolic  curve  fit  of  the 
experimental  data,  an  integration  routine  using  the  coefficients 
A,  B  and  C  of  the  curve  fit  was  derived.  The  integracion 
procedure  was  designed  to  Integrate  any  integrand  that  could  be 
approximated  by  the  curve  fit.  In  the  case  of  the  experimental 
data,  the  Integrand  was  written  as 

Cp(s)  -  As2  +  Bs  +  C  .  (A. 1.3) 

The  integral  I(s)  was  evaluated  as  the  following  definite  integral, 

8  2  8  2 

I(s)  -  f  (As2  +  Bs  +  C)ds  -  [y  S3  +  |  S2  +  Cs]  I  .  (A. 1.4) 

81  81 

Therefore, 

I(s)  -  (s 2  -  sjf-  (s2  +  +  S2)  +|  (s2  +  Sj)  +  C] 

(A. 1.5) 

which  represents  the  value  of  the  integral  I(s)  across  the 
specified  time  interval. 

When  computing  the  bubble-time  parameter  from  Eq.  (2.8)  the 
integrand,  one  writes  f(s)  as 


Substitution  of  the  curve  fit  for  Cp(s)  into  Eq.  (A. 1.6)  produced 

tabulated  values  of  the  integrand  f(s)  for  different  arc  length 

positions  along  the  body.  Application  of  the  parabolic  curve  fit 

routine  to  Eq.  (A. 1.6)  produces  a  piecewise  parabolic  approximation 

of  the  integrand  f(s)  which  is  then  integrated  exactly  as  C  (s) 

P 

was  integrated  using  Eq.  (A. 1.5). 

Because  the  code  is  required  to  produce  an  output  point  for 
every  input  point  and  to  handle  an  even  or  odd  number  of  data 
points,  two  separate  Integration  formulae  are  used.  The  first 
integration  formula,  1^,  Integrates  across  the  first  two  points 
in  a  triad  of  data  points  while  the  second  formula,  I2,  integrates 
across  the  entire  triad  of  data  points.  By  integrating  across 
consecutive  triads  of  data,  i.e., 

[x(  l),x(2),x(3) ]  ,  (x(3) ,x(4) ,x(5) ] ,  ...  , 

one  can  get  a  total  sum  for  the  integral  found  in  Eq.  (2.8)  which 
allows  one  to  find  the  bubble  time  for  different  arc  length 
positions.  Figure  A1  shows  how  the  integral  is  computed  for  an 
even  or  odd  number  of  data  points.  When  n  *  5,  there  is  an  even 
number  of  intervals  and  1^  and  1^  are  computed  across  two  triads 
of  data  points.  When  n  *  6,  there  is  an  odd  number  of  intervals 
and  1^  arad  are  computed  across  two  triads  of  data  points  and 
across  one  group  of  two  data  points.  Since  requires  three  data 


points  across  which  It  Is  calculated  we  must  set  at  x(7)  equal 
to  zero.  This  procedure  Is  used  across  the  entire  set  of  data  to 
produce  the  required  1^  for  1  ■  l,2,...,n.  Figure  A2  is  the  flow 
chart  for  the  integration  procedure.  It  satisfies  all  requirements 
put  forth  In  Chapter  II  but  it  does  not  describe  the  procedures 
for  computing  the  coefficients  A,  B  and  C.  The  specifics  of  these 
calculations  are  described  in  Section  3  of  Appendix  A. 

A2.  Evaluation  of  the  Stagnation  Point  Singularity 

Flow  about  a  symmetrical  hemispherical  headform  produces  a 
stagnation  point  at  the  nose  of  the  headform  characterized  by  a 
zero  fluid  velocity  and  a  maximum  pressure.  The  history  of  a  very 
small  flaccid  air  bubble  traveling  close  to  the  stagnation 
streamline  would  show  a  constantly  decreasing  radius  up  to  the 
stagnation  point  where  a  minimum  radius  occurs.  Past  th 
stagnation  point  along  a  meridian  on  the  body  the  flaccid  nubble 
would  experience  an  increasing  radius.  It  should  be  noted  that 
we  speak  of  the  bubble  as  traveling  close  to  the  stagnation 
streamline  and  not  exactly  on  the  stagnation  streamline  because 
it  would  take  an  infinite  amount  of  time  to  reach  the  stagnation 


point  as  it  travels  from  infinity.  Also,  along  the  body  arc  the 
flow  velocity  is  zero,  but  the  bubble  is  traveling  in  the 
boundary  layer  and  not  on  the  body. 
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Figure  A2.  Flow  Chart  Representing  the  Algorithm  by  Which  the 
Integral  is  Calculated  for  an  Even  or  Uneven  Number 
of  Data  Points. 


If  we  define  a  pressure  coefficient,  ,  as 


C 

P 


P 


P 

o 


l/2pVo2 
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(A. 2.1) 


we  could  expect  to  see  a  minimum  bubble  radius  at  a  maximum  pressure 

coefficient  of  unity  and  a  maximum  bubble  radius  at  a  minimum 

pressure  coefficient,  C  ,  provided  that  no  cavitation  takes 

^min 

place.  The  definition  of  the  pressure  coefficient  leads  to  a 
singularity  when  computing  the  dimensionless  translation  time  for 
a  bubble  on  the  stagnation  streamline.  As  determined  earlier,  the 
dimensionless  time  was  written  as 

T1  .  i  -  1.2 . .  (A. 2. 2) 

0  V  „r3  /I  -  Cp(B) 

where  the  singularity  in  the  integrand  is  an  integrable  square  root 
singularity.  Since  the  computer  code  was  not  designed  to  integrate 
across  a  singularity,  an  analytic  technique  was  used  on  the 
interval  containing  the  singularity  while  the  remainder  of  the  data 
was  Integrated  by  use  of  the  computer  code. 

The  interval  containing  the  singularity  was  composed  of  the 
following  points: 


£l 

1.0000 

0.9895 

0.9791 


s_ 

0.00000 

0.00279 

0.00557 


Since  the  computer  code  produced  a  piecewise  parabolic  approximation 


of  the  experimental  data,  we  will  write  the  equation  for  Cp  across 
the  interval  containing  the  singularity  as 

Cp(s)  -  As2  +  Bs  +  1  .  (A. 2. 3) 


Evaluation  of  the  coefficients  across  the  interval  containing  the 
singularity  gives 


A  -  -  0.87307 
B  -  -  3.74416  . 

Substituting  Cp(s)  into  the  integral  one  gets 


which  when  solved  using  math  tables  [12)  one  finds  I(s)  »  0.000983. 
This  value  of  I  is  then  added  onto  the  values  of  the  integral  for 
successive  intervals  as  computed  by  the  computer  code.  Figure  A3 
and  Fig.  A4  show  the  integrand,  f(s). 


f  (a) 


1 

✓  l  -  C  (s) 
P 


and  the  integral,  I(s), 


I(s)  -  /  -  'ds 

s ,  /I  -  C  (s  ) 

1  P 

each  plotted  against  the  body  arc  length  s.  It  is  clearly  seen 
that  the  singularity  at  the  nose  of  the  body  causes  the  integrand 
to  increase  to  infinity  but  that  the  integral  at  the  singularity 
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A3.  Derivation  of  the  Computer  Code  Formulae 

Having  tabulated  information  on  the  pressure  distribution 
around  a  hemispherical  headform,  one  must  be  able  to  integrate  the 
data  based  on  the  bubble  time  integral  discussed  in  Chapter  2. 

It  was  stated  before  that  these  data  are  unevenly  spaced  and  may 
be  odd  or  even  in  number.  Also,  since  we  desire  one  output  point 
for  each  input  point,  we  need  two  integration  fornulae  to  do  this. 
For  tabulated  data  of  the  form  (x^,f(x^))  where  i  -  1,2,. . .  ,n  is 
being  the  total  number  of  data  points,  the  integration  equation  will 
be  developed  for  a  triad  of  data  points  and  then  used  across  each 
successive  triad  until  all  data  have  been  used. 

The  procedure  for  defining  the  integration  formulae  begins  by 
calculating  a  parabolic  curve  fit  for  each  triad  of  data.  The 
coefficients  A,  B  and  C  are  calculated  using  Eqs.  (A. 1.1)  and 
(A. 1.2)  as  discussed  in  Appendix  A.  By  expanding  these  equations 
for  i  ■  1,2,3  one  gets 


p„(.) 


3 

E 

i-1 


3 

x  (s  -  S  ) 

j-1  J 

ifj _ 

n 

x  (s  -  s  ) 

j-1  3 

i*j 


,  ,  (s  “  s2Hs  ‘  *3) 

f(V  (.J  -s2)(3l  -.3) 


w  ^  (s  -  ^H3  -  s3)  w  ^  (s  -  ^H3  -  s2) 

+  f  S2  (s2  -  SjJ(s2  -  s3J  +  f(S3'  (s3  -  s1j(s3  -  s2j 

Factoring  like  powers  of  s  one  gets  an  equation  resembling  a 
parabola  which  looks  like 
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The  parabolic  form  of  curve  fit  can  now  be  integrated  across  the 
appropriate  limits  to  obtain  the  two  integration  formula  required 
for  obtaining  an  output  point  for  each  input  point  when  the  number 
of  data  points  can  be  odd  or  even  and  the  data  are  unevenly 
spaced*  Therefore,  we  define  the  two  integration  formulae  as 
Ij  and  I ^  where 

S2 

I,  *  /  p  (s  )ds 
l  n 

8i 

and 

S3 

I-  -  /  P  (s)ds  . 

L  n 


Upon  evaluation  of  the  definite  integrals,  one  can  write 


These  integration  formulae  can  be  used  for  any  data  which  are 
tabulated  in  the  form  of  s^,  f(s^)  and  can  be  modeled  by 
successive  parabolic  curve  fits  according  to  the  parameters 
A,  B,  C  and  A  defined  earlier  in  this  section.  One  must  keep 
in  mind  that  the  integration  routine  is  only  as  good  as  the 


ESSSSGSl  v&sxs&tti-  m&maa 


curve  fit  since  the  coefficients  of  the  curve  fit  play  a  major  role 
in  the  result  of  the  integration.  The  remainder  of  Appendix  A 
will  deal  primarily  with  the  accuracy  of  the  curve  fit  and 
integration  formulae. 

A4.  Accuracy  of  the  Parabolic  Curve  Fit 

When  using  any  type  of  numerical  approximation,  it  is  important 
that  the  user  be  aware  of  the  limitations  and  accuracy  of  the 
particular  approximation  technique  being  used.  Failure  to  use  an 
approximation  technique  with  suitable  accuracy  could  result  in 
erroneous  calculations  and  analysis. 

In  this  work,  it  was  imperative  that  the  approximation 
technique  used  to  calculate  the  piecewise  parabolic  curve  fit 
and  numerical  integration  were  extremely  accurate.  The  curve 
fit  approximation  was  compared  against  a  sample  function  and  the 
accuracy  was  checked  for  different  interval  lengths. 

The  sample  function  used  was  f(x)  *  sin  x.  The  first 
interval  length  used  was  Dx  ■  0.157  for  x  ■  0  to  x  ■  t.  The 
second  interval  length  to  be  used  was  Dx  ■  0.314  for  x  -  0  to 
x  ■  ir  while  a  third  interval  length  of  Dx  ■  0.628  was  used 
across  x  *  0  to  x  *  i. 

Execution  of  the  code  produced  the  coefficients  of  the 
parabolic  curve  fit  which  was  compared  to  the  sample  function. 

When  using  the  shortest  interval,  Dx  *  0.157,  the  code 
produced  a  curve  fit  accurate  to  five  decimal  places.  Doubling 


the  interval  length  produced  a  curve  fit  accurate  to  four 


decimal  places.  And  when  the  interval  length  was  made  four  times 
the  first  interval  or  Dx  -  0.628,  the  accuracy  dropped  to  three 
decimal  places. 


Due  to  the  irregularities  of  the  data,  several  points  were 
interpolated  between  the  original  data  points.  This  effectively 


reduced  the  interval  lengths  so  that  the  curve  fit  was  accurate 


to  five  decimal  places.  Therefore,  the  method  described  in  the 
previous  sections  of  the  Appendix  produced  a  curve  fit  which 
allowed  an  accurate  fit  of  the  data.  Section  1  of  Appendix  A 
addressed  the  problem  im  more  detail  but  the  accuracy  of  the 
results  will  be  given.  Again,  using  a  sample  function  to  test  the 
integration  procedure  against  the  parabolic  coefficients  of  the 
curve  fit  needed  to  be  calculated.  For  the  same  function 


f(x)  ■  sin  x  from  x  *  0  to  x  •  it  the  value  of  the  definite  integral 
is  exactly  equal  to  2.  Using  the  parabolic  coefficients  in  the 
integration  procedure  gave  the  value  of  the  integral  equal  to 
1.999977.  This  calculation  was  performed  with  an  interval  length 
of  Dx  ■  0.157.  Further  calculations  were  performed  on  even  and 
odd  numbers  of  data  points  with  longer  and  shorter  interval  and 
it  was  found  that  the  accuracy  of  the  integration  procedure 
paralleled  the  accuracy  of  the  curve  fit.  Thus,  when  using 
the  curve  fit  routine  and  integration  procedure  in  this  work, 
the  interval  lengths  were  short  enough  to  produce  five  decimal 


place  accuracy 


APPENDIX  B.  DATA  TABULATION 


BL .  Conversion  of  the  Axial  Length  Data  to  Arc  Length  Data 

Based  on  the  geometry  of  the  headform  set  forth  In  Fig.  1, 
the  following  equations  are  used  to  convert  the  axial  length, 
x/D,  which  Is  the  parameter  that  Holl  and  Carroll  [1]  measured  the 
0^  data  against,  to  a  dimensionless  arc  length  s.  The  parameter 
s  was  used  throughout  this  work  to  describe  the  dimensionless  arc 
length  along  the  surface  arc  of  the  body. 

For  0°  <  9  <  90° 

s  -  y  cos  1  (l  -  2(x/D) ) 

For  9  >  90° 

s  -  0.785  +  | 

It  should  be  noted  that  for  0  >  90°  the  dimensionless  arc  length 
s  Increases  the  same  as  the  parameter  x/D  due  to  the  presence 
of  the  cylindrical  afterbody.  Therefore,  the  value  of,  x/D  is 
simply  added  on  to  the  value  of  s  beyond  0  ■  90°.  Section  2 
of  Appendix  B  contains  a  tabulation  of  s  versus  x/D  with  the 
corresponding  averaged  values  of  C  . 
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B2.  Cp  versus  8  (averaged) 


x/D 

s 

S’ 

0.000 

0.000 

1.0000 

0.280 

0.557 

-  0.6168 

0.335 

0.617 

-  0.7343 

0.390 

0.674 

-  0.7812 

0.430 

0.715 

-  0.7432 

0.465 

0.748 

-  0.6578 

0.480 

0.765 

-  0.6370 

0.500 

0.785 

-  0.6163 

0.515 

0.800 

-  0.6128 

0.530 

0.815 

-  0.5518 

0.545 

0.830 

-  0.4183 

0.560 

0.845 

-  0.3303 

0.575 

0.860 

-  0.2988 

0.625 

0.910 

-  0.2307 

0.675 

0.960 

-  0.2073 

0.725 

1.010 

-  0.1805 

0.775 

1.060 

-  0.1620 

0.825 

1.110 

-  0.1443 

